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ABSTRACT 


To  recover  3-D  structure  from  a  shaded  and  textural  surface  image  involving 
textures,  neither  the  Shape-from-shading  nor  the  Shape-from-texture  analysis  is 
enough,  because  both  radiance  and  texture  information  coexist  within  the  scene 
surface.  A  new  3-D  texture  model  is  developed  by  considering  the  scene  image  as  the 
superposition  of  a  smooth  shaded  image  and  a  random  texture  image.  To  describe  the 
random  part,  the  orthographical  projection  is  adapted  to  take  care  of  the  non-isotropic 
distribution  function  of  the  intensity  due  to  the  slant  and  tilt  of  a  3-D  texture  surface, 
and  the  Fractional  Differencing  Periodic  (FDP)  model  is  chosen  to  describe  the 
random  texture,  because  this  model  is  able  to  simultaneously  represent  the  coarseness 
and  the  pattern  of  the  3-D  texture  surface,  and  enough  flexible  to  synthesize  both 
long-term  and  short-term  correlation  structures  of  random  texture.  Since  the  object  is 
described  by  the  model  involving  several  free  parameters  and  the  values  of  these 
parameters  are  determined  directly  from  its  projected  image,  it  is  possible  to  extract 
3-D  information  and  texture  pattern  directly  from  the  image  without  any  pre¬ 
processing.  Thus,  the  cumulative  error  obtained  from  each  pre-processing  can  be 
minimized.  For  estimating  the  parameters,  a  hybrid  method  which  uses  both  the  least 
square  and  the  maximum  likelihood  estimates  is  applied  and  the  estimation  of 
parameters  and  the  synthesis  are  done  in  frequency  domain.  Among  the  texture 
pattern  features  which  can  be  obtained  from  a  single  surface  image,  Fractal  scaling 
parameter  plays  a  major  role  for  classifying  and/or  segmenting  the  different  texture 


patterns  tilted  and  slanted  due  to  the  3-dimensional  rotation,  because  of  its  rotational 
and  scaling  invariant  properties.  Also,  since  the  Fractal  scaling  factor  represents  the 
coarseness  of  the  surface,  each  texture  pattern  has  its  own  Fractal  scale  value,  and 
particularly  at  the  boundary  between  the  different  textures,  it  has  relatively  higher 
value  to  the  one  within  a  same  texture.  Based  on  these  facts,  a  new  classification 
method  and  a  segmentation  scheme  for  the  3-D  rotated  texture  patterns  are  developed. 
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CHAPTER  1 

INTRODUCTION  AND  OVERVIEW 


1.1.  Introduction 

An  important  task  in  computer  vision  is  the  recovery  of  3-D  scene  information 
from  single  2-D  images.  However,  this  is  an  ill-posed  problem  [13, 178],  because  we 
have  only  2-D  information  but  try  to  extract  3-D  information  from  it.  Therefore,  we 
need  additional  information  for  this  missing  dimension  to  recover  3-D  scene 
information.  To  date,  researchers  have  suggested  two  different  ways  to  handle  this 
lack  of  information.  The  first  one  is  the  use  of  shading  information  in  the  image 
[30,31,73,75, 100, 144, 172],  Shading  information  tells  us  the  direction  of  the  light 
source  and  the  surface  orientation  of  the  object  surface.  Thus,  by  formulating  the 
reflectance  map  function  which  shows  scene  radiance  as  a  function  of  the  surface 
gradient  and  the  distribution  of  light  sources,  we  can  extract  3-D  surface  information 
from  image  data.  The  second  method  is  the  use  of  texture  information 
1 6, 82,  109, 135,207,223],  Since  texture  gradients  behave  like  intensity  gradients,  the 
shape  of  a  surface  can  be  inferred  from  the  pattern  of  a  texture  on  the  surface  by 
applying  statistical  texture  analysis. 

However,  for  describing  a  natural  scene  image,  each  of  the  above  approaches 
have  their  own  limitations.  Shading  information  is  meaningful  only  under  the 
assumption  that  the  surface  is  smooth  enough  to  have  clear  radiance  information,  even 
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though  this  situation  is  rarely  encountered  in  practice.  Thus,  instead  of  having  clear 
radiance  information,  the  scene  image  has  texture  pattern  information  due  to  the 
complexity  of  the  surface  in  most  cases.  This  encourages  us  to  study  the  texture 
pattern  of  a  3-D  surface.  Therefore,  the  modeling  and  analysis  of  texture  patterns 
which  contain  3-D  information  of  the  surfaces  will  be  the  main  focus  of  this  thesis. 

3-D  texture  analysis  involves  at  least  the  several  following  difficulties.  First,  for 
the  recovery  of  3-D  shape  from  the  natural  scene,  a  texture  model  itself  will  not  be 
enough  to  represent  the  whole  surface,  because  texture  analysis  requires  the  surface  to 
be  relatively  complex  to  have  a  statistical  property.  However,  some  parts  of  the 
surface  may  not  contain  clear  information  about  the  texture  pattern  due  to  dark 
shading  or  bright  radiance.  In  other  words,  we  may  have  a  situation  such  that  the 
radiance  information  is  more  dominant  than  the  texture  information.  Therefore,  our 
model  should  have  an  ability  to  handle  both  the  texture  and  radiance  information  at 
the  same  time.  Second,  since  the  observed  texture  image  is  actually  the  image 
obtained  by  projecting  the  surface  image  to  the  image  plane,  this  projected  random 
texture  pattern  does  not  have  the  stationary  property  any  more,  even  though  the 
original  surface  normal  image  can  be  represented  by  a  stationary  random  texture 
pattern.  Thus,  the  conventional  2-D  texture  models  under  the  assumption  of 
stationarity  can  not  be  used.  Third,  for  the  classification  and  texture  segmentation 
purpose,  our  model  should  have  a  3-D  rotational  invariant  property,  since  the 
observed  texture  pattern  may  have  various  looks  depending  on  the  viewer’s  direction. 
Thus,  the  classification  and/or  segmentation  scheme  for  the  3-D  rotated  texture  pattern 
should  have  the  flexibility  to  treat  different  looking  textures  as  being  the  same  without 
losing  the  accuracy  of  the  processes. 

In  this  thesis,  to  solve  the  difficulties  discussed  above,  a  composite  model  of 
Shape-from-shading  and  Shape-from-texture  is  developed  to  represent  a  3-D  surface 
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image  considering  the  scene  image  as  the  superposition  of  a  smooth  shaded  image  and 
a  random  texture  image.  The  orthographical  projection  is  adapted  to  take  care  of  the 
non-isotropic  distribution  function  due  to  the  slant  and  tilt  of  a  3-D  texture  surface, 
and  the  Fractional  Differencing  Periodic  model  is  chosen  because  this  model  is  able  to 
simultaneously  represent  the  coarseness  and  the  pattern  of  the  3-D  texture  surface 
with  the  fractional  differencing  parameters  c,  d  and  the  frequency  parameters  g)j,  0)2  • 
For  the  classification  and  segmentation  purpose,  these  fractional  differencing 
parameters  play  an  important  role,  because  those  parameters  are  known  to  be 
rotational  and  scaling  invariant.  Thus,  combining  these  parameters  and  directional 
frequency  parameters  in  the  fractional  differencing  periodic  function,  we  can  have  the 
flexibility  to  handle  the  rotated  and  projected  texture  pattern  and  the  accuracy  of  the 
classification.  For  estimating  the  parameters,  a  hybrid  method  which  uses  both  the 
least  square  and  the  maximum  likelihood  estimates  is  applied  and  the  estimation  and 
the  synthesis  are  done  in  frequency  domain. 

1.2.  Modeling  Of  A  Surface  Image 

Modeling  of  a  3-D  surface  image  can  be  broken  down  into  two  main  categories, 
Shape-ffom-shading  and  Shape-from-texture.  The  Shape-from-shading  model  uses 
the  reflectance  map  which  shows  scene  radiance  as  a  function  of  the  surface  gradient 
and  the  distribution  of  light  sources  to  extract  3-D  surface  information  from  image 
data  [30,99].  On  the  other  hand,  the  Shape-from-texture  model  uses  the  texture 
pattern  instead  of  shading  to  extract  3-D  structure  [82,109,135].  Since  texture 
gradients  behave  like  intensity  gradients,  the  shape  of  a  surface  can  be  inferred  from 
the  pattern  of  a  texture  on  the  surface  by  applying  statistical  texture  analysis. 
However,  for  describing  a  natural  scene  image,  each  of  the  above  approaches  has  its 
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own  limitation.  The  Shape-from-shading  model  is  applicable  only  under  the 
assumption  that  the  surface  is  smooth  enough  to  have  clear  radiance  information, 
while  the  Shape-from-texture  model  requires  the  surface  to  be  relatively  complex  so 
that  texture  information  can  be  extracted.  Thus,  neither  model  is  suitable  to  represent 
a  natural  scene,  because  both  radiance  and  texture  information  coexist  within  the 
surface  of  a  natural  scene.  Therefore,  a  robust  technique  is  needed  to  handle  this 
shortcoming.  Recently,  the  fractal  scaling  parameter  was  introduced  to  measure  the 
coarseness  of  the  surface,  and  applied  to  represent  the  natural  scene  surface 
[47, 132, 151, 173].  However,  this  fractal  model  is  also  not  enough  to  represent  the 
real  3-D  texture  image,  because  even  though  two  surfaces  are  estimated  to  have  the 
same  fractal  scales,  these  surfaces  can  have  different  texture  patterns.  Therefore,  in 
this  thesis,  the  Fractional  Differencing  Periodic  model  is  chosen  to  represent  the  3-D 
surface  image,  because  this  model  is  able  to  simultaneously  represent  the  coarseness 
and  the  pattern  of  the  3-D  texture  surface  with  the  fractional  differencing  parameters 
c,  d  and  the  frequency  parameters  (Oi ,  o^- 

1.3.  Fractional  Differencing  Model 

As  mentioned  previously,  the  fractional  differencing  model  has  an  ability  to 
simultaneously  represent  the  coarseness  [175]  and  the  pattern  of  the  3-D  texture 
surface  with  the  fractional  differencing  parameters  c,  d  and  the  frequency  parameters 
coi ,  (1)2.  Also,  it  has  the  property  of  being  flexible  enough  to  synthesize  both  long-term 
and  short-term  correlation  structures  of  random  texture  depending  on  the  values  of  the 
fractional  differencing  parameters  [89, 101].  For  estimating  the  parameters, 
comparing  with  the  fractional  Brownian  random  process  model  [151],  the  fractional 
differencing  model  has  a  simple  estimation  scheme  sharing  the  same  properties  of  the 
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fractional  Brownian  process,  because  while  the  fractional  Brownian  process  is  a 
continuous  process  which  follows  a  certain  probability  distribution,  the  fractional 
differencing  model  is  a  discrete  process  which  has  a  linear  function  of  parameters.  A 
hybrid  method  which  uses  both  the  least  square  and  the  maximum  likelihood 
estimates  is  applied  and  the  estimation  and  the  synthesis  are  done  in  the  frequency 
domain  [671. 

1.4.  Shape  From  A  Shaded  And  Textured  Surface  Image 

In  this  thesis,  a  composite  model  of  Shape-from-shading  and  Shape-from-texture 
is  developed  to  represent  a  3-D  surface  image  considering  the  scene  image  as  the 
superposition  of  a  smooth  shaded  image  and  a  random  texture  image,  that  is,  the 
deterministic  function  x(l\,li)  and  the  random  function  y{l j,^).  Then,  the 
orthographical  projection  is  adapted  to  take  care  of  the  non-isotropic  distribution 
function  due  to  the  slant  and  tilt  of  a  3-D  texture  surface.  Thus,  the  Orthographically 
Projected  Fractional  Differencing  Periodic  (OPFDP)  model  is  chosen  because  this 
model  is  able  to  simultaneously  represent  the  coarseness  and  the  3-D  rotated  pattern  of 
the  surface  with  the  fractional  differencing  parameters  c,  d,  the  frequency  parameters 
C0j ,  o>2,  and  the  relationship  between  different  directions  of  3-D  surface.  Since  the 
object  is  described  by  a  model  involving  several  free  parameters  and  the  values  of 
these  parameters  are  determined  directly  from  its  projected  image,  it  is  possible  to 
extract  3-D  information  and  texture  pattern  directly  from  the  given  intensity  values  of 
the  image  without  any  pre-processing.  Thus,  the  cumulative  error  obtained  from  each 
pre-processing  can  be  minimized.  For  estimating  the  parameters,  a  hybrid  method 
which  uses  both  the  least  square  and  the  maximum  likelihood  estimates  is  applied  and 
the  estimation  and  the  synthesis  are  done  in  frequency  domain  based  on  the  local 
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patch  analysis.  By  using  this  model,  the  integrability  problem  wnich  might  occur  in 
spatial  domain  analysis  can  be  avoided,  because  only  one  inverse  Fourier  transform 
needs  to  be  taken  at  the  end  of  the  procedure  to  get  the  whole  image. 

1.5.  Classification  Of  3-D  Rotated  Textures 

The  classification  problem  can  be  stated  as  an  allocation  of  observed  texture 
image  data  to  the  one  of  the  pre-defined  texture  classes.  These  texture  classes  can  be 
described  by  texture  features,  and  then  texture  features  can  be  the  parameters  in 
stochastic  model  [46,53,92],  or  structural  model  [142,216].  Thus,  the  key  step  in  the 
classification  process  is  the  choice  of  a  set  of  features  which  can  reduce  the  dimension 
of  the  image  data  to  a  computationally  reasonable  amount  of  data.  The  features  should 
be  simple  and  easy  to  extract  from  the  given  data  while  preserving  the  classifying 
information  present  in  the  data. 

Most  classification  schemes  which  have  been  suggested  to  date  are  under  the 
assumption  that  the  test  sample  data  possesses  the  same  surface  orientation  as  the 
training  sample  data.  Thus,  if  the  orientation  of  test  image  is  different  from  the 
training  sample  data,  for  example,  in  case  of  a  rotated  image,  the  classification 
performs  poorly.  This  reduces  the  flexibility  of  those  classification  schemes. 
However,  most  natural  texture  images  which  we  can  encounter  in  practice  are 
representing  the  texture  on  the  3-D  surface,  thus,  the  observed  image  is  a  projected 
surface  image  onto  the  2-D  image  plane  with  a  2-D  rotated,  or  3-D  slanted  and  tilted 
texture  pattern.  Therefore,  sometimes  it  is  desirable  for  the  classification  scheme  to 
have  the  flexibility  that  it  can  classify  even  rotated  or  scaled  texture  to  the  original 
class  of  it.  This  is  a  good  indication  why  it  is  so  important  to  have  the  rotational  and 
scaling  invariant  features  in  our  model.  In  this  thesis,  a  classification  method  which 
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can  handle  arbitrary  3-D  rotated  samples  of  textures  is  developed,  i.e.,  the  accuracy  of 
classification  is  not  affected  by  the  3-D  rotation  of  the  test  texture.  This  classification 
scheme  is  based  on  a  two-level  hierarchical  structure.  In  the  first  level,  a  3-D 
rotational  invariant  feature,  fractal  scale  c,  is  extracted  from  the  first-order  fractional 
differencing  model  by  applying  a  least-square  estimation  method,  and  this  feature  is 
used  to  classify  the  test  texture  image  to  a  class  whose  members  are  sharing  a  similar 
value  of  fractal  scale.  And  in  the  second  level,  the  members  of  the  class  are  further 
classified  to  the  final  desired  subclasses  with  other  texture  pattern  features,  coj  and  0)2, 
which  are  extracted  from  the  second-order  projected  fractional  differencing  model  by 
applying  a  hybrid  method  of  the  least-square  estimation  and  the  maximum  likelihood 
estimation.  As  a  result,  this  multi-level  classification  scheme  saves  a  reasonable 
amount  of  processing  time  without  losing  the  accuracy  of  the  classification. 

1.6.  Organization  Of  The  Thesis 

Various  applications  of  the  Fractional  Differencing  model  have  been  investigated 
to  represent  3-D  texture  pattern  through  this  thesis.  An  important  aim  of  this  study  is 
to  develop  the  mathematical  model  suitable  for  the  3-D  surface  image  which  the 
radiance  and  texture  information  coexist  in.  The  Orthographically  Projected  Fractional 
Differencing  (OPFD)  model  developed  here  performs  very  well  to  represent  the 
texture  pattern  on  the  3-D  surface,  because  of  the  rotational  and  scaling  invariant 
parameters  in  it.  This  rotational  and  scaling  invariant  property  of  these  parameters  has 
been  successfully  applied  to  segment  or  classify  the  rotated  and  slanted  texture  plane 
in  the  rest  of  the  chapters. 
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The  organization  of  the  thesis  is  as  follows.  In  chapter  2,  two  different  categories 
of  3-D  surface  model,  Shape-from-shading  and  Shape-from-texture,  are  discussed,  and 
several  typical  methods  of  each  categories  are  compared  with  their  simulation  results. 
Two  projection  methods,  orthographical  and  perspective  projection,  are  also  presented 
to  represent  the  distortion  due  to  the  tilt  and  the  slant  of  the  surface.  In  chapter  3,  the 
Fractional  Differencing  model  suggested  in  chapter  2  to  represent  the  3-D  texture  is 
discussed  in  detail,  and  its  estimation  scheme  based  on  the  Least-Square  and 
Maximum  likelihood  estimation  methods  is  presented.  Chapter  4  presents  a  composite 
model  of  Shape-from-shading  and  Shape-from-texture  to  extract  the  3-D  structure 
from  the  surface  image  which  contains  the  radiance  information  and  the  texture 
information  at  the  same  time.  This  suggested  model  is  directly  applied  to  the  given 
image  without  any  pre-processing,  and  as  a  result  of  this,  the  errors  which  might  result 
from  each  pre-processing  are  not  cumulated.  In  chapter  5,  a  classification  scheme  of 
the  3-D  rotated  textures  is  developed  based  on  the  fractal  scale.  This  fractal  scale  is 
known  to  be  a  rotational  and  scaling  invariant  parameter,  and  can  be  extracted  by 
fitting  the  given  tilted  or  slanted  texture  image  to  the  proposed  Fractional  Differencing 
model.  A  multi-level  structure  of  the  classification  structure  is  also  introduced  to 
reduce  the  processing  time.  Finally  in  chapter  7,  the  conclusion  of  this  study  and  the 
suggested  future  research  are  presented. 


9 


CHAPTER  2 

MODELING  OF  A  SURFACE  IMAGE 


2.1.  Introduction 

Modeling  of  a  3-D  surface  image  can  be  broken  down  into  two  main  categories, 
Shape-ffom-shading  and  Shape-from-texture.  The  Shape-from-shading  model  uses 
the  reflectance  map  which  shows  scene  radiance  as  a  function  of  the  surface  gradient 
and  the  distribution  of  light  sources  to  extract  3-D  surface  information  from  image 
data  [30,73,75, 100, 144, 172].  On  the  other  hand,  Shape-from-texture  model  uses  the 
texture  pattern  instead  of  shading  to  extract  3-D  structure  [6, 109, 135,207,216,221]. 
Since  texture  gradients  behave  like  intensity  gradients,  the  shape  of  a  surface  can  be 
inferred  from  the  pattern  of  a  texture  on  the  surface  by  applying  statistical  texture 
analysis.  However,  for  describing  a  natural  scene  image,  each  of  the  above 
approaches  has  its  own  limitation.  The  Shape-from-shading  model  is  applicable  only 
under  the  assumption  that  the  surface  is  smooth  enough  to  have  clear  radiance 
information,  while  the  Shape-from-texture  model  requires  the  surface  to  be  relatively 
complex  so  that  texture  information  can  be  extracted.  Thus,  neither  model  is  suitable 
to  represent  a  natural  scene,  because  both  radiance  and  texture  information  coexist 
within  the  surface  of  a  natural  scene.  Therefore,  a  robust  technique  is  needed  to 
handle  this  shortcoming.  Recently,  the  fractal  scaling  parameter  was  introduced  to 
measure  the  coarseness  of  the  surface,  and  applied  to  represent  the  natural  scene 
surface  [47,132,151,173].  However,  this  fractal  model  is  also  not  enough  to 
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represent  the  real  3-D  texture  image,  because  even  though  two  surfaces  are  estimated 
to  have  the  same  fractal  scales,  these  surfaces  can  have  different  texture  patterns. 
Therefore,  in  this  thesis,  the  Fractional  Differencing  Periodic  model  is  chosen  to 
represent  the  3-D  surface  image,  because  this  model  is  able  to  simultaneously 
represent  the  coarseness  and  the  pattern  of  the  3-D  texture  surface  with  the  fractional 
differencing  parameters  c,  d  and  the  frequency  parameters  coj ,  ©2. 

2.2.  Shape-From-Shading 


Pioneering  work  on  the  inference  of  shape-from-shading  was  done  by  Horn  and 
his  co-workers  [30,99, 100],  To  extract  3-D  shape  from  a  single  2-D  image,  they  used 
the  reflectance  map,  which  shows  the  intensity  of  the  image  as  a  function  of  the 
surface  gradient  and  the  illumination  direction.  (Figure  2.1  shows  the  relationship 
between  the  different  directions  and  angles.) 

pcos xL  sin  OL+qs\ n xL  si  n  Gj, +cos  G/, 


/(/ 1,/2)  = 


(P2V+ 1) 


1/2 


kR{p,q) 


(2.2.1. a) 


=  sincjcosxcosx^sino^  +  sinosinxsinx^sina^  +  cosacosa^ 


(2.2. 1  .b) 


where  p  =  -^j-H(li,l2),  q  = 

/?(p,q) :  Reflectance  map  function 

:  3-D  shape  function  from  the  viewing  direction, 
x,  a  :  Tilt,  slant  of  the  surface 

'■  Tilt,  slant  of  the  illumination  direction 

Here,  the  relationship  between  x,  a  and  p,  q  are 


p  -  tanacosx,  q  =  tanasint 


(2.2.2. b) 


Figure  2.2-a,b  were  simulated  from  the  equation  (2.2.1),  assuming  that  the  slope 
values  p,  q  in  directions  of  / 1 , 1 2  and  the  illumination  direction  were  given. 


(a)  (b) 

Figure  2.2:  Intensity  images  of  a  sphere  from  the  various  directions  of  light 
source,  (a)  =  0.0,T^  =  0.0  (b)  =  -0.66, Xi  =  -0.66  [rad],  (Images  are 

simulated  by  (2.1.2.1a)  and  intensity  values  are  normalized  between  0  and 
255.) 


Construction  of  3-D  shape  can  be  achieved  by  solving  p,  q  in  terms  of  /(/  j  ,(2)  at 
each  point,  and  integrating  those  values.  However,  the  final  integrated  shape  can  be 
different  from  the  original  shape,  due  to  the  cumulation  of  estimation  errors.  To  avoid 
this  type  of  error,  Horn  and  Brooks  [100]  developed  a  calculus  of  variation  method  to 
estimate  the  surface  orientation  values,  and  Pentland  [172]  suggested  local  shape 
analysis  which  deals  with  the  only  local  areas  instead  of  a  whole  image.  However, 
Pentland’s  technique  has  severe  trouble  in  integrating  all  local  areas.  Recently, 
Pentland  [169[  developed  another  technique  to  solve  the  integrability  problem.  He 
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suggested  analysis  in  the  frequency  domain,  instead  of  in  the  spatial  domain.  By  using 
this  method,  the  integrability  problem  can  be  avoided,  because  only  one  inverse 
Fourier  transform  needs  to  be  taken  at  the  end  of  procedure.  However,  since  the 
calculation  of  convolution  is  required  in  the  frequency  domain  to  handle  simple 
multiplication  operation  of  spatial  domain,  calculation  will  be  complicated.  Horn  and 
Brooks’  method  was  improved  by  Frankot  and  Chellappa  [76]  recently,  by  adding  one 
more  constraint  of  the  integrability  in  addition  to  the  smoothness  constraint.  The 
estimation  error  and  the  integrability  problem  can  be  handled  by  the  image  model 
suggested  in  this  thesis,  because  our  model  has  a  constraint  which  is  from  the  texture 
pattern  and  is  analyzed  in  frequency  domain. 

2.2.1.  Estimation  Of  Illumination  Direction 

The  first  step  in  estimating  local  surface  orientation  is  the  determination  of  the 
illumination  direction  for  the  surface,  L,  and  the  constant  ^.(i  for  a  particular 
estimation  neighborhood.  Estimation  methods  of  illumination  direction  L  were 
suggested  by  Pentland  [171],  and  C-H  Lee  [144]  in  different  ways. 

2.2.1. 1.  Pentland’s  Method 

Pentland’s  approach  was  based  under  the  assumption  that  the  surface  normal  of 
each  local  patch  is  isotropically  distributed  within  a  scene.  Assuming  that  the 
distribution  of  the  surface  normal  is  known  and  the  intensity  value  is  measurable  for 
different  directions  within  the  image,  one  could  estimate  the  illumination  direction  L 
using  a  least  square  estimation  procedure  [171],  The  solution  of  this  approach  follow. 
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Theorem  2.1:  Let  di(l  i,l2)  be  the  first  derivative  of  image  intensity  at  a  particular 
point  in  the  direction  ( dl\ ,  dl2),  and  dij  be  the  average  value  of  di  over  the  j-th  patch 
in  direction  ( dl \j,dl2j).  Then, 


dx  1 
dx  2 

dxn 


dl\\  dl 21 
dl  12  dl 22 


dl  in  dl  2/1 


LL 

4 


The  proof  is  in  Appendix  A. 


Theorem  2.2  :  The  illumination  directions  for  1 1 ,  /  2, 1 3  direction  are 


Li  I - ; 


where  k=  J(di2)-(di ) 

di 2  :  the  average  value  of  di2  for  the  whole  image. 
di :  the  average  value  of  di  for  the  whole  image, 
and 


-1  ^2  -1 

~  tan  l(— ),  aL  =cos  lLh 


where  x 1  :  Tilt  angle  of  illumination  direction 

<3l  :  Slant  angle  of  illumination  direction 
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The  proof  is  in  Appendix  B. 

Pentland’s  method  is  a  very  good  tilt  estimator,  as  well  as  a  very  good  slant 
estimator  if  the  slant  is  below  40  degrees.  However,  due  to  the  small  variance  of 
d i(lj  ,12),  this  method  results  in  a  large  error  in  the  slant  estimation  whenever  the  slant 
is  above  50  degrees.  This  experimental  result  has  been  confirmed  by  Lee  [144]  and 
Ferrie  &  Levine  [73] 

2.2. 1.2.  Lee’s  Method 

Lee  suggested  a  different  statistical  approach  to  estimate  the  illumination 
direction  [144],  By  assuming  that  the  slant  a  is  uniformly  distributed  over  the  sphere 
surface,  the  probability  density  function  for  the  slant  G  can  be  obtained  by  (sin  a)/2  k. 
To  determine  the  expected  value  of  image  brightness  i(\\  ,12),  we  should  integrate  only 
over  the  illuminated  region,  because  not  all  of  the  sphere  is  illuminated.  Thus,  the 
region  of  integration  should  be  set  based  on  the  foreshortening  factor.  Also,  this 
region  of  integration  can  be  different  under  the  condition  when  the  self-shadowed 
area",  is  included  in  the  computation  of  the  average.  If  we  include  the  self-shadowed 
areas  in  the  computation,  we  must  divide  the  integral  by  the  whole  area,  n,  of  the  disc 
that  is  the  projection  of  the  hemisphere.  If,  on  the  other  hand  we  do  not  include  the 
self-shadowed  areas,  we  divide  by  the  area  (;t/2)(1+cosGl)  of  the  projection  of  the 
illuminated  part  of  the  hemisphere.  Thus,  from  the  integration,  we  can  get  the 
expected  values  of  i(li  J2)  and  / 2 (1 1 ,12)  as  follows: 


Theorem  2.3  :  The  slant  of  illumination  direction  Gl  satisfies 
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Thus,  the  slant  Cl  can  be  estimated  from  above  equations  by  taking  the  average 
value  of  intensities  over  the  whole  image  for  E[/]  and  E[/2].  The  tilt  of  illumination 
direction  Tl  can  be  also  estimated  from  Theorem  4.1  of  [144]. 

EI-gpiGt.12)] 

TL  =  tan-1( - ^ - )  (2.2. 1.2.1) 

E[— i(li,l2)3 

Here,  the  estimated  values  E(-)  are  taken  over  the  whole  image  by  taking  the  average 
value  of  the  function. 

Lee’s  method  has  several  advantages  over  Pentland’s  method.  First,  calculation 
is  much  simpler,  because  calculations  are  required  only  once  over  the  whole  image. 
Second,  we  can  calculate  the  slant  of  illumination  direction  directly  from  two 
equations  without  knowing  the  value  of  fyt.  Experimental  results  of  Lee’s  method  are 
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known  to  be  superior  to  Pentland’s  for  the  estimation  of  both  tilt  and  slant.  This  has 
been  confirmed  by  Lee  [144]  and  Ferrie  &  Levine  [73].  Because  of  these  advantages, 
Lee’s  method  was  used  to  determine  the  illumination  direction  for  the  construction  of 
3-D  shape  in  this  thesis. 

2.2.2.  Estimation  Of  Surface  Orientation 

The  extraction  of  3-D  surface  orientation  from  single  2-D  images  is  an  ill-posed 
problem,  because  we  have  only  2-D  information  but  try  to  extract  3-D  information 
from  it.  Therefore,  we  need  additional  information  for  this  missing  dimension  to 
recover  3-D  scene  information.  To  date,  several  approaches  have  been  suggested  by 
researchers  to  give  an  additional  constraint,  such  as  the  integrability,  smoothness,  etc. 
The  following  sections  will  discuss  several  approaches  among  them  in  detail. 

2.2.2.I.  Horn  And  Brooks’s  Method 

Horn  and  Brooks  [30, 100]  developed  a  calculus  of  variation  method  to  estimate 
surface  orientation.  In  the  presence  of  noise,  the  real  values  may  not  be  the  same  as 
the  estimated  values  that  satisfy  the  image  irradiance  equation  exactly.  There  will, 
however,  be  a  surface  that  minimizes  the  integral  of  the  square  of  the  error  between 
the  expected  values  and  the  real  values.  Thus,  the  search  for  a  function  that  minimizes 
an  integral  of  this  error  was  taken  to  be  the  major  concern  of  this  calculus  of 
variations.  However,  this  problem  is  an  ill-posed  problem,  because  there  are  typically 
an  infinite  number  of  surface  satisfying  this  equation.  Therefore,  the  equation  needs  an 
additional  constraint  to  have  a  unique  solution.  Horn  and  Brooks  proposed  the 
additional  constraint  of  a  smoothness  criterion  [30].  Surface  orientation  can  be 
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obtained  by  minimizing  the  following  cost  function 

jjm  ,\2)-R(p.q))2+H(^YH(k  ,l2))2+2(^^tf(li  ,12))2 

+  (-ii-//(l1,l2))2Ml1Jl2  (2.2.2. 1.1) 

ol2 

Horn’s  method  has  good  estimation  results  in  the  presence  of  noise.  However,  when 
the  surface  is  not  relatively  smooth  enough,  the  result  from  this  method  can  possess 
the  integrability  problem.  This  lack  of  integrability  can  be  found  from  the  simulation 
examples  in  Section  2.2.3. 

2.2.2.2.  Chellappa  And  Frankot’s  Method 

Frankot  and  Chellappa’s  method  [76,78]  is  based  on  the  calculus  of  variation 
method  of  Horn  and  Brooks.  However,  this  method  deals  with  one  more  constraint 
about  the  integrability.  As  discussed  before,  Horn  &  Brooks  added  the  smoothness 
constraint  to  yield  the  unique  solution,  but  that  method  can  have  the  lack  of 
integrability.  Thus,  Frankot  and  Chellappa  proposed  an  additional  constraint  to 
enforce  the  integrability.  This  integrability  can  be  achieved  by  satisfying  the  following 
requirement. 

ai5irw<1-w=3i5irw<,"l2)  <2-2'2'21) 

There  are  many  conceivable  ways  of  enforcing  above  equation.  One  of  them  is  the 
minimization  the  following  distance  measure. 

JJ((p-E[p])2+(q-E[ql)2)dMl2  (2.2.2.2.2) 

Let  C  be  the  coefficients  of  the  Fourier  series  expansion  of  H(\\ ,12),  Q,  and  be  the 
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Fourier  coefficients  for  the  expectation  values  of  p,  q  (slopes  in  directions  of  lj,^). 
Then,  the  above  cost  function  can  be  minimized  by  taking 

C(co)  = - 5 5 -  (.3) 

<V+0V 

Thus,  the  construction  of  3-D  shape  can  be  achieved  by  simply  taking  the  inverse 
Fourier  transform  on  C(co). 

Frankot  and  Chellappa’s  method  is  a  good  estimator  of  surface  orientation. 
From  the  smoothness  constraint  and  enforcing  integrability,  the  general  shape  of  the 
3-D  object  can  be  obtained.  However,  due  to  too  much  smoothing,  the  detail 
information  about  the  surface  may  be  lost. 

2.2.2.3.  Pentland’s  Method 


The  Shape-from-shading  problem  is  known  to  be  mathematically  equivalent  to  a 
nonlinear  first-order  partial  differential  equation  in  surface  elevation.  Therefore,  it  is 
very  difficult  to  obtain  the  closed  form  solution.  For  this  reason,  Pentland  [169] 
proposed  a  linear  approximation  method  to  estimate  the  values,  using  a  Taylor  series 
expansion  in  p,  q  variables  up  to  the  second  order.  Thus, 


z  (1 1 ,12)  =  cosOL+pcos'CLSinaL+qsin'TLsinoL- 


COSCTl 


(p2+q2)  (2.2.2.3.1) 


The  corresponding  DFT  of  /(li  ,12)  i-s  as  follows,  (after  deleting  the  constant  term 
cosgl) 

cosGl  ^ 

fi(k1,k2)  =  costLsinGLfp(k1,k2)+sinxLsinGLfq(k1,k2) - - — (fp  ®fP+fq  ®fq) 
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(2.2.23.2) 

Then,  fp(k!,k2)  and  fq(k!,k2)  can  be  represented  as  the  function  of  f//(kj ,k2),  by  the 
approximation  of  the  first  derivative  of  H{ li,l2)  function.  Thus,  f4(ki,k2)  can  be 
represented  as  a  function  of  f//(k!,k2)  only,  and  the  height  function  H{ li,l2)  can  be 
constructed  by  simply  taking  the  inverse  Fourier  transform  of  the  intensity  function 

/O1.I2). 

f  P(ki  ,k2)  =  7sin(27t-^-)f  //  (kx  ,k2)  (2.2.2.33) 

Similarly, 

fq(k„k2)  =  ysin(27t~)f x-(k,  ,k2)  (2.2.23.4) 

Linder  the  condition  |p|,  |q|  <  1,  the  linear  term  dominates.  Thus, 
f  ,(kj  ,k2)  =  cosxLsincLf  p(k,  ,k2)  +  sintLsinaLf  q(k]  ,k2) 

k]  ^2 

=  7'sinaL(sin(27t— )costL  +  sin(27t— )sinTL)f  //(kj  ,k2)  (2.2.23.5) 

Pentland’s  method  has  a  nice  mathematical  formulation  because  in  frequency  domain, 
the  equation  has  only  one  unknown  variable.  Thus,  this  equation  has  a  unique 
solution.  However,  the  first-order  linear  approximation  yields  a  big  estimation  error, 
while  the  second-order  approximation  requires  convolution  in  frequency  domain  and 
requires  too  much  computation.  This  results  will  be  confirmed  from  the  examples  in 
Section  2.23. 
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2.2.3.  Simulation  Results 

A  128  x  128  real  image  of  a  part  of  a  tree  in  a  grass  meadow  is  considered 
(Figure  2.3).  From  this  image,  the  estimation  values  of  the  illumination  direction  <Tl, 

A  A 

tl  was  determined  to  be  <Jl=  -1.3384,  Tl=  -4.257 899e-02,  by  applying  Lee’s  method. 
Following  Figure  2.4-a,b  show  the  results  of  the  height  functions  constructed  from 
Horn  and  Brooks’s  and  Frankot  and  Chellappa’s  algorithms,  respectively.  From  these, 
we  can  see  that  Frankot’s  algorithm  gives  a  nice  and  smooth  height  function,  while 
Horn’s  gives  a  rough  surface  but  detailed  information  of  the  surface.  The  experiments 
were  repeated  for  the  smoothed  part  images  which  were  obtained  by  different  sizes  of 
smoothing  windows.  These  experiments  could  give  the  idea  of  how  the  smoothing 
affects  the  results.  Figure  2.5-a,b,c  depict  the  pan  images  which  are  smoothed  by  th? 
different  sizes  of  window,  and  Figure  2.6-a,b,  Figure  2.7-a,b,  and  Figure  2.8-a,b  also 
depict  the  corresponding  height  functions  from  each  cases  by  Horn’s  and  Frankot’s, 
respectively.  From  these  experiments,  we  can  see  that  for  the  image  smoothed  by  the 
relatively  large  size  of  window,  Horn’s  and  Frankot’s  are  getting  closer.  Pentland’s 
algorithm  was  not  considered  here,  because  of  the  poor  performance  from  the  linear 
approximation.  Figure  2.9  shows  the  resulting  height  function  constructed  from  the 
low-frequency  linear  approximation  method  of  Pentland  for  Figure  2.5-b. 
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Figure  2.3  :  A  digitalized  tree  image  sized  by  512  x  512 


(a) 


Figure  2.4:  The  height  function  of  Figure  2.3:  constructed  by  (a)  Flom’s  method, 
(b)  Frankot’s  method:  Surface  orientation  for  each  pixel  was  obtained  from 
minimizing  the  cost  function  (2.2.2. 1.1)  and  (2. 2.2. 2. 2),  respectively. 
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(b) 


Figure  2.4,  continued 
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(b) 


(c) 


Figure  2.5:  A  part  of  the  tree  image  (Figure  2.3):  (a)  original,  (b)  Image  obtained 
after  convoluting  a  5  x  5  smoothing  window  to  (a),  (c)  Image  obtained  after 
convoluting  a  9x  9  smoothing  window  to  (a). 
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(a) 


Figure  2.6:  The  height  function  of  Figure  2.5-a:  constructed  by  (a)  Horn’s 
method,  (b)  Frankot’s  method:  Surface  orientation  for  each  pixel  was  obtained 
from  minimizing  the  cost  function  (2. 2. 2. 1.1)  and  (2. 2.2.2. 2),  respectively. 


(a) 


Figure  2.7:  The  height  function  of  Figure  2.5-b:  constructed  by  (a)  Horn’: 
method,  (b)  Frankot’s  method:  Surface  orientation  for  each  pixel  was  obtaima 
from  minimizing  the  cost  function  (2.2.2. 1 . 1 )  and  (2. 2. 2. 2. 2;,  respectively. 


Figure  2.8:  The  height  function  of  Figure  2.5-c:  constructed  by  (a)  Horn's 
method,  (b)  Frankot’s  method:  Surface  orientation  for  each  pixel  was  obtained 
from  minimizing  the  cost  function  (2.2.2  l.H  and  (2.2  2  2  2),  respectively 


2.3.  Shape-From-Texture 


The  Shape-from-texture  problem  has  been  studied  extensively  to  recover  the  3-D 
information  from  a  two  dimensional  image  [6, 109, 135, 136, 156,221]  since  Gibson 
first  proposed  the  texture  density  gradient  as  the  primary  basis  of  surface  perception 
by  humans  in  1950.  Motivated  to  simulate  human  perception,  he  confines  the  surfaces 
under  consideration  to  only  those  whose  orientations  could  be  easily  perceived  from 
the  texture  by  human  eye.  However,  if  the  true  texture  is  not  isotropic  and  has  a 
preferred  orientation,  it  mimics  a  projected  image  and  thus  makes  it  difficult  to  detect 
the  true  orientation.  Therefore,  to  formalize  the  shape-from-texture  problem  requires  a 
projection  model  for  the  image  formation  system.  Up  to  now,  two  kinds  of  projections 
called  ‘orthographical  projection’  and  ‘perspective  projection’  have  been  used  in  most 
of  cases.  The  orthographical  projection  can  be  obtained  by  projecting  the  object 
surface  to  the  image  plane  parallelly  without  considering  the  average  distance  from 
the  camera,  while  in  the  perspective  projection  model,  the  scene  depth  is  relative  to 
the  distance  between  the  object  and  the  camera.  Thus,  the  perspective  transformation 
yields  orthographical  projection  as  a  special  case  when  the  viewpoint  is  at  a  point 
infinitely  far  from  the  object  surface.  Another  model  required  to  formalize  the  shape- 
from-texture  is  a  random  field  model.  This  model  should  represent  the  statistical 
property  of  the  texture  on  the  surface  properly,  and  its  parameter  estimation  scheme 
should  exist.  Various  types  of  random  field  models  have  been  proposed  to  date  to 
represent  the  surface  covered  by  texture.  Among  them,  the  AR(Auto-Regressive) 
model  [41, 122]  is  known  to  be  simple  enough  to  estimate  and  synthesize  the  texture 
plane  whose  statistical  pattern  is  isotropically  distributed  over  the  plane.  Another 
noticeable  model  for  the  3-D  textured  surface  is  the  Fractal  model  [133,173].  The 
fractal  model  contains  a  fractal  scaling  parameter.  This  fractal  scaling  parameter  can 
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be  a  measure  of  the  roughness  of  the  surface  and  can  be  considered  as  a  scale  and 
rotational  invariant  parameter.  However,  these  models  have  their  own  limitations  in 
representing  the  true  texture  surface.  The  AR  model  does  not  have  a  long-term 
persistent  memory  property,  thus,  it  is  not  proper  for  the  texture  pattern  which 
contains  a  long-term  memory  in  a  certain  direction,  such  as  a  tree-bark  image.  On  the 
other  hand,  the  Fractal  model  has  only  three  variable  parameters,  mean,  variance  and 
fractal  scaling.  Those  are  not  flexible  enough  to  model  the  wide  range  of  situations 
encounted  in  practice.  This  conflict  can  be  solved  in  the  fractional  differencing 
periodic  model  which  is  suggested  in  this  paper.  That  is,  me  tractional  differencing 
periodic  model  contains  four  parameters,  two  frequency  parameters,  which  are  similar 
to  the  texture  pattern  parameters  in  AR  model,  and  two  other  fractional  differencing 
parameters,  which  are  corresponding  to  the  fractal  scaling  parameter  in  the  fractal 
model.  Therefore,  the  fractional  differencing  periodic  model  gives  more  flexibility  of 
modeling. 

2.3.1.  Projection 

As  mentioned  in  previous  section,  the  projection  model  is  basic  and  necessary  to 
formalize  the  shape-from-texture  problem  to  represent  the  slanted  and  tilted  texture 
surface.  To  date,  there  are  two  different  projection  models  which  have  been  used  in 
most  cases.  Those  are  ‘Orthographical  projection’  and  ‘Perspective  projection’.  Detail 
discussion  on  these  will  be  given  in  the  following  sections. 
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2.3.1. 1.  Orthographical  Projection 


Consider  a  plane  with  a  texture  on  it.  and  take  the  nipn^  coordinate  system.  Put 
a  line  passing  through  the  origin,  and  let  T  be  the  angle  made  from  the  tri)  -axis.  Rotate 
the  plane  around  the  line  by  angle  a  and  project  the  rotated  plane  orthographically 
onto  the  original  plane  (Figure  2.10). 


m)  orl2 


Figure  2.10:  Coordinate  transformation  of  the  orthographic  projection  (2.3. 1.1.4) 


Thus,  a  new  coordinate  system  from  the  viewing  direction,  lj-lj,  can  be  obtained 
from  the  following  two  coordinate  transformations. 


Vi' 

cost  -sinT 

*'2 

sinx  cost 

m2 

(2.3.1. 1.1) 


and 
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’ll' 

1  0 

V,' 

>2 

0  coso 

1'2 

(2.3.1. 1.2) 


Hence,  as  is  well  known,  the  coordinate  transformation  of  the  orthographic 
projection  between  the  mi -m2  system  and  lpb  system  can  be  given  as  follows. 


* 

“ 

r 

“ 

“ 

-1 

' 

h 

cosx  -sinx 

1 

0 

cosx  -sinx 

m. 

12 

sinx  cosx 

0 

COSO 

j 

sinx  cosx 

m2 

(2.3.1.1.3-a) 


Thus, 


■>  •  2 
cos  x+cososin  x 

(l-coso)sinxcosx 


(l-coso)sinxcosx 

sin2x+cosocos2x 


(2.3.1.1.3-b) 


m, 

1 

sin2x+cosocos2x 

(coso-1  jsinxcosx 

r 

ll 

mi 

COSO 

(coso-1  jsinxcosx 

cos2x+cososin2x 

b 

One  grid  pattern  image  (Figure  2.11 -a)  was  considered  to  demonstrate  this 
orthogonal  projection.  The  coordinate  transformation  was  taken  to  this  image  with 

O  =  — -,  x  =  —  and  a=  — ,  x  —  — — .  Figure  2.1 1  -b  and  Figure  2.1 1  -c  depict  these 

8  4  8  4  b  h 

transformation  respectively. 


(a) 


(b) 


(c) 


Figure  2.11:  2-D  grid  pattern  images:  (a)  Original,  (b)  Image  obtained  after 
projecting  the  image  (Figure  2.1 1 -a)  orthographically,  with  o  =  n/4  and  x  =  rc/8  in 
(2.3.1. 1.4),  (c)  Image  obtained  after  projecting  the  image  (Figure  2.11  -a) 
orthographically,  with  o  =  rc/4  and  x  =  — tx/8  in  (2.3. 1.1. 4). 
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2.3. 1.2.  Perspective  Projection 


The  perspective  projection  acts  like  a  pinhole  camera  in  that  the  image  results 
from  projecting  scene  points  through  a  single  point  onto  an  image  plane. 


Figure  2.12  :  Perspective  projection 


As  in  Figure  2.12,  if  the  point  of  projection  corresponds  to  a  viewpoint  behind 
the  image  plane  and  the  image  occurs  right  side  up,  the  viewpoint  is  +f  on  the  z  axis, 
with  z  =  0  plane  being  the  image  plane  upon  which  the  image  is  projected.  Thus,  as 
the  image  object  approaches  the  viewpoint,  its  projection  gets  bigger.  In  Figure  2.12, 
y',  the  projected  height  of  the  object,  is  related  to  its  real  height  y,  its  position  z,  and 
the  focal  length  f  by 

y'  =  (f/f— z)y  (2.3. 1.2.1) 


Similarly 
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x'  =  (f/f-z)x 


(2.3. 1.2.2) 


The  projected  image  has  z  =  0  everywhere.  However,  projecting  away  the  z 
component  is  best  considered  a  separate  transformation;  the  projective  transform  is 
usually  thought  to  distort  the  z  component  just  as  it  does  the  x  and  y.  Perspective 
distortion  thus  maps  (x,  y,  z)  to 


(x',  y\  z')  =  ( 


fx  fy 


fz 


f-  z  ’  f-  z  ’  f-  z 


) 


(2.3. 1.2.3) 


Notice  that  the  perspective  transformation  yields  orthographical  projection  as  a  special 
case  when  the  viewpoint  is  the  point  at  infinity  in  the  z  direction.  Then  all  objects  are 
projected  onto  the  viewing  plane  with  no  distortion  of  their  x  and  y  coordinates. 
Another  point  about  the  perspective  transformation  is  that  if  z  is  a  known  function  of 
x,  y  and  the  surface  orientation  a  and  %,  the  result  of  the  transformations  in  (2.3. 1.2.1) 
and  (2.3. 1.2.2)  is  a  nonhomogeneous  planar  texture.  This  is  important  because  under 
orthographical  projection  a  homogeneous  texture  on  a  slanted  and  tilted  plane  was 
transformed  to  another  homogeneous  texture  due  to  the  linearity  of  the  transformation 
(2. 3. 1.1. 4),  while  under  perspective  projection  it  is  not  a  homogeneous  texture  any 
more.  Thus,  the  analysis  in  the  frequency  domain  by  the  discrete  Fourier  transform  is 
no  longer  valid. 

To  overcome  this  difficulty  associated  with  the  perspective  projection,  we  can 
approximate  the  perspective  projection  by  an  affine  transformation  which  is  suggested 
by  Aloimonos  and  Swain  [6].  The  approximation  is  done  by  dividing  the  projection 
process  into  two  steps.  The  first  step  is  projecting  the  local  plane  Q  with  orientation 
given  by  oj  onto  the  plane  T  which  is  parallel  to  image  plane  I.  This  projection  is 
performed  parallel  to  the  ray  OG,  where  G  is  the  center  of  local  plane.  The  second 
step  is  projecting  this  plane  T  perspectively  onto  the  image  plane  I.  Since  the  plane  T 
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is  parallel  to  the  image  plane,  this  perspective  transformation  is  just  a  reduction  by  a 
factor  of  1/(3.  Figure  2.13  depicts  this  relationship. 

Thus,  the  transformation  from  surface  plane  coordinate  (irq  ,nri2)  to  image  plane 
coordinate  (I1T2)  with  the  two  step  projection  process  is  given  by  the  following  affine 
transformation  [6]. 


-1+pA  q(p-l-A) 

’ll' 

1 

Vi+p2  V(i+p2)(i+p2+q2) 

mi 

h 

”  P 

pB  qB-p2-l 

m2 

_Vl+p2  V(l+p2)(l+p2-t-q2)  _ 

Here,  the  point  (A,B,-1)  is  the  mass  center  of  the  image  texel,  and  p,  q  are  the  same  as 
defined  in  (2.2.1).  Therefore,  the  non-homogeneous  texture  due  to  perspective 
projection  can  be  approximated  by  piecewise  homogeneous  ones,  and  each  local 
homogeneous  texture  patch  is  easily  synthesized  by  using  the  above  linear  affine 
transformation. 

2.3.2.  Random  Texture  Analysis 

For  the  surface  which  is  covered  by  a  texture  pattern,  or  is  relatively  complex, 
the  random  field  model  can  be  applied  over  the  surface  normal  plane  for 
■'.pproximating  the  surface  image.  However,  differently  from  the  regular  2-D  texture 
analysis,  3-D  textural  surface  image  analysis  can  be  enhanced  by  considering  the 
fractional  differencing  parameter,  (which  is  the  ‘fractal  scaling  parameter’  in  the 
terminology  of  Pentland  [173, 176])  which  indicates  the  roughness  of  the  surface,  in 
addition  to  the  texture  pattern.  In  other  words,  a  3-D  textured  surface  can  be 
represented  by  the  fractal  scaling  and  the  texture  pattern  parameters  in  a  certain 
random  field  model.  Therefore,  the  model  based  on  2-D  texture  pattern  only,  such  as 


Figure  2.13  :  Approximation  method  of  perspective  projection 


the  Auto-Regressive(AR)  model  [122],  or  the  model  based  on  2-D  fractional 
Brownian  motion  process  with  one  fractal  scaling  parameter,  such  as  Pentland’s 
model  [176],  has  limited  modeling  ability.  That  is  because  the  AR  model  or  the  fractal 
model  can  represent  only  the  texture  pattern  or  the  roughness  of  the  surface,  thus, 
either  model  is  not  flexible  enough  to  represent  the  wide  range  of  different  texture 
patterns  encountered  in  practice,  especially  non-stationary  random  texture.  This 
limitation  can  be  overcome  by  using  the  2-D  FDP  (Fractional  Differencing  Periodic) 
model  which  contains  both  texture  pattern  and  fractal  scaling  parameters  [67]. 

2.3.2.I.  At*lo-Regressive(AR)  Model 

The  autoregressive  model  is  one  of  the  classical  short  correlation  models.  This 
model  is  simple  but  has  good  modeling  performance.  Assume  that  the  image  intensity 
follows  the  autoregressive  model.  Let  (I1J2)  be  an  index  for  the  coordinate  location 
and  i(li , I2)  be  the  intensity  at  the  coordinate  Then  the  auto-regressive  model  is 

represented  by  the  following  equation. 

i(l1,l2)  =  eTz(l1,l2)  +  C(li,l2)  (2.3.2.1.1) 

where  0  is  parameter  vector  and  z(lj  ,12)  is  a  vector  which  consists  of  intensities  of  the 
neighbor  pixels  and  unity.  {  C0i>l2)  }  is  a  two  dimensional  white  noise  sequence 
with  variance  p. 

For  example,  in  case  of  the  Causal  Auto-Regressive  (CAR)  model  as  a  special 
case  of  the  AR  model,  the  parameter  vector  0  is  4  x  1  column  vector  and  zO^L)  is 
another  4  x  1  column  vector  which  contains  three  causal  neighbors  and  unity  as  the 
elements. 
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In  3-D  textural  surface  image,  the  fractional  differencing  parameter,  which  is 
‘fractal  scaling  parameter’  in  the  terminology  of  Pentland  [173,176],  indicates  the 
roughness  of  the  surface.  That  is,  as  the  value  of  the  fractal  scaling  parameter 
increases,  the  model  represents  a  3-D  surface  textured  more  roughly.  Typical  2-D 
fractional  Brownian  motion  surface  iflj,^),  which  was  suggested  by  Pentland, 
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satisfies  the  following  relationship.  Define  Ai^r  =  |  i(li2»l22)  -  KI11.I21)  I  > 

Ar  =  V0i2-lii)2  +  (l22-l2i)2-  Then 

E(AiAr)-ArF  (2. 3. 2.2.1) 

The  fractal  dimension  of  the  fractional  Brownian  motion  surface  is  then  3  -  F 
[151,46].  However,  this  model  is  not  proper  to  present  3-D  texture  pattern  because 
the  fractal  of  surface  is  assumed  to  be  spatially  isotropic,  thus,  the  fractional  Brownian 
motion  process  with  one  fractal  parameter  is  not  flexible  enough  to  represent  the  wide 
range  of  different  texture  patterns  encountered  in  practice,  especially  non-stationary 
random  texture.  Also,  comparing  with  other  model-based  model  such  as  AR  model, 
this  model  is  not  practical  to  be  used  because  it  is  difficult  to  estimate  the  fractal 
scaling  parameter  from  the  given  intensity  image  directly. 

2.3.2.3.  Fractional  Differencing  Model 

The  1-D  Fractional  Differencing  model  was  suggested  by  Hosking  [101], 
generalizing  the  well-known  ARIMA  model  of  Box  &  Jenkins  [23],  which  was 
originally  designed  to  model  a  non-stationary  random  process.  This  1-D  model  was 
extended  to  2-D  case  by  Kashyap  and  Lapsa  [116]  later.  A  typical  second-order 
Fractional  Differencing  Periodic  (FDP)  model  is  as  follows. 

C 

iOi  .I2) =  (l-2cosa>!  Z!-l+  Z!-2)  2 

_  d_ 

•  (1-2cosco2  z2_1  +z2-2)  2  C(li,l2)  (2.3.2.3.1) 


for  1 1 , 12  =0,1, . N-l  . 
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The  corresponding  DFT  of  this  function  is 

-1  kl  A  kl  C 

-j2n-n-  “J 4nTT  -T 

I(k!,k2)  =  (l-2coscoie  N+e  N)  2 

-3271^  -j47i-i 

•(l-2coso)2e  N+e  N)2W(k1,k2),  (2.3.23.2) 

where  z;  is  the  delay  operator  associated  with  1;,  £(li,l2)  is  an  i.i.d.  Gaussian 
sequence,  and  W(k}  ,k2)  is  the  corresponding  DFT. 

This  model  has  four  different  parameters,  c  and  d  for  the  fractal  scales  and  CO] 
and  (1)2  for  the  frequencies  of  pattern  in  the  direction  of  lj  and  12,  respectively.  Thus, 
this  model  represents  the  roughness  of  the  surface  and  the  pattern  of  the  texture  image 
at  the  same  time  even  with  the  different  values  for  the  direction  of  lj ,  12  separately. 
Also  notice  that  this  random  process  model  has  the  property  of  being  flexible  enough 
to  explain  both  the  long-term  and  short-term  correlation  structure  of  a  time  series 
depending  on  the  values  of  the  fractional  differencing  parameter  c,  d,  and  it  shares  the 
basic  properties  with  Fractional  Brownian  motion  defined  by  Mandelbrot  [151],  More 
detailed  discussion  about  estimation  scheme  and  synthesis  of  this  model  will  be  given 
in  the  next  section. 

2.3.3.  Simulation  Results 

For  the  experiment,  one  tree  bark  texture  on  the  surface  of  a  tree  is  considered 
(Figure  2.14-a).  As  seen  from  Figure  2.14-a,  the  pattern  of  tree  bark  texture  has 
different  values  of  frequency  and  roughness  in  the  horizontal  and  vertical  directions. 
Thus,  it  cannot  be  represented  by  either  Pentland’s  model  or  a  2nd-order  AR  model, 
because  Pentland’s  model  can  have  only  various  fractal  scales  and  a  2nd-order  AR 


46 


model  has  only  directional  frequencies.  Figure  2.14-b.c  show  the  limitation  of 
Pentland’s  Fractal  models  with  fractal  scale  c=0.9  and  c=1.5,  respectively,  and  Figure 
2.14-d  shows  the  limitation  of  a  2nd-order  AR  model  with  the  directional  frequencies 
co i=0.25  and  0)2=0.025  in  direction  of  mj,  m2. 

On  the  other  hand,  Figure  2.15  shows  how  well  the  Fractional  Differencing 
Periodic  (FDP)  model  fits  to  the  tree  bark  texture.  Here,  the  directional  frequencies 
coi,  0)2  and  the  fractal  scales  c,  d  in  FDP  model  are  chosen  as  0.25,  0.025,  1.5,  0.8, 
respectively. 

2.4.  Conclusions 

Various  models  to  represent  3-D  surface  have  been  considered.  The  Shape- 
from-shading  model  is  based  on  the  reflectance  map  function  which  shows  scene 
radiance  as  a  function  of  the  surface  gradient  and  the  distribution  of  light  sources  to 
extract  3-D  surface  information.  The  Shape-from-texture  model  uses  the  texture 
pattern  instead  of  shading  to  extract  3-D  structure.  Then,  the  fractal  model  was 
discussed  considering  a  fractal  scale  to  represent  the  roughness  of  the  surface.  Thus, 
several  estimation  schemes  for  the  direction  of  the  light,  the  projection  methods,  and 
2-D  statistical  texture  models  have  been  considered.  However,  for  describing  a  natural 
scene  image,  the  above  approaches  have  their  own  limitations.  The  Shape-from- 
shading  model  is  applicable  only  under  the  assumption  that  the  surface  is  smooth 
enough  to  have  clear  radiance  information,  the  Shape-from-texture  model  requires  the 
surface  to  be  relatively  complex  so  that  texture  information  can  be  extracted,  and  the 
fractal  model  with  one  fractal  scaling  parameter  does  not  have  enough  flexibility  to 
represent  the  various  texture  patterns  which  we  can  encounter  in  practice.  To  give 
more  flexibility  to  represent  a  real  surface  image  and  to  have  an  ability  to  handle  the 


47 


1 5 


Figure  2.13:  Synthesized  surface  shapes  over  64x64  sized  normal  plane  patch  from 
(a)  A  tree  image  whose  surface  is  covered  by  the  tree  bark  texture,  (b)  Pentland’s 
model  with  the  fractal  scale  c=0.9,  (c)  Pentland’s  model  with  the  fractal  scale 
c=1.5,  (d)  2nd-order  AR  model  with  the  frequency  CO]  =0.25,  g>2=0.()25  in 
directions  of  m( ,  mi,  respectively. 


smooth  surface  and  the  textured  surface  simultaneously,  the  Fractional  Differencing 
model  was  suggested  in  this  chapter. 


(a) 


(b) 


Figure  2.14:  Surface  shapes  obtained  from  Fractional  Differencing  Periodic  model 
(2.3.2.3.1):  (a)  Synthesized  surface  shape  over  64x64  sized  normal  plane  patch 
with  the  frequency  w^O.25,  0)2=0.025  and  the  fractal  scale  c=1.5,  d=0.8  in 
direction  of  rrq ,  m2,  (b)  corresponding  surface  image. 
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CHAPTER  3 

FRACTIONAL  DIFFERENCING  MODEL 
AND  ESTIMATION 


3.1.  Introduction 

As  mentioned  previously,  the  fractional  differencing  model  has  an  ability  to 
simultaneously  represent  the  coarseness  and  the  pattern  of  the  3-D  texture  surface 
with  the  fractional  differencing  parameters  c,  d  and  the  frequency  parameters  CD!,  (D2- 
Also  it  has  the  property  of  being  flexible  enough  to  synthesize  both  long-term  and 
short-term  correlation  structures  of  random  texture  depending  on  the  values  of  the 
fractional  differencing  parameters  [  101 J.  For  estimating  the  parameters,  comparing 
with  the  fractional  Brownian  random  process  model  [151],  the  fractional  differencing 
model  has  a  simple  estimation  scheme  sharing  the  same  properties  together,  because 
while  the  fractional  Brownian  process  is  a  continuous  process  which  follows  a  certain 
probability  distribution,  the  fractional  differencing  model  is  a  discrete  process  which 
has  a  linear  function  of  parameters.  A  hybrid  method  which  uses  both  the  least  square 
and  the  maximum  likelihood  estimates  is  applied  and  the  estimation  and  the  synthesis 
are  done  in  frequency  domain  [67,  127], 
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3.2.  Fractional  Brownian  Model 

Brownian  motion  is  a  continuous  time  stochastic  process  B  ( t )  with  independent 
Gaussian  increments  and  spectral  density  oT2.  Its  derivative  is  the  continuous-time 
white  noise  process,  which  has  constant  spectral  density.  Fractional  Brownian  motion, 
Bp(t),  introduced  by  Mandelbrot  and  Van  Ness  [  1 5 1 J,  is  a  generalization  of  these 
process  as  follows.  Let  0  <  F  <  1  and  let  b q  be  an  arbitrary  real  number.  Define 
Bf(  0)  =  b0, 

,  o 

BF{t)  -  Bf(  0)  =  —  {  f  [(r-5)f-1/2  -  (-sf-m]dB  (s) 

r (F  + 1/2) 

i 

+  \{t-s)F~mdB  (s) },  for  t  >  0,  (3.2.1) 

o 

and  similarly  for  t  <  0.  Notice  that  for  F  =  1/2,  this  definition  becomes  that  of 
classical  Brownian  motion  B  (r).  This  Fractional  Brownian  motion  has  the  following 
basic  properties  1101], 

Property  1:  Fractional  Brownian  motion  with  parameter  F,  usually  0  <  F  <  1,  is  the 
(y  -  F)th  fractional  derivative  of  Brownian  motion. 

Property  2:  The  spectral  density  of  fractional  Brownian  motion  is  proportional  to 
(0_2/r_1. 

Property  3:  The  covariance  function  of  fractional  Brownian  motion  is  proportional 
to  | k  \2F~2. 

Thus,  the  derivative  of  fractional  Brownian  motion,  B'p(t),  may  also  be  thought  of  as 
the  (y-F)th  fractional  derivative  of  continuous-time  white  noise,  to  which  it  reduces 
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when  F  =  -.  One  of  the  possible  discrete-time  analogues  of  this  continuous-time 

fractional  noise  is  discrete-time  fractional  Gaussian  noise,  which  is  defined  to  be  a 
process  whose  correlation  function  is  the  same  as  that  of  the  process  of  unit 
increments  ABp(t)  =  Bp(t)  -Bp(t- 1)  of  fractional  Brownian  motion. 

However,  the  modeling  ability  of  fractional  Brownian  motion  noise  has 
sometimes  been  claimed  to  be  inferior  to  that  of  other  processes  such  as  the  high-order 
Moving  Average  (MA)  process.  This  is  because  of  that  fractional  Brownian  model  has 
only  three  variable  parameters,  mean,  variance  and  F,  and  those  are  not  flexible 
enough  to  model  the  wide  range  of  randomness  encounted  in  practice.  Also,  it  is 
relatively  difficult  to  estimate  parameter  F  and  to  synthesize  with  it,  because  of  the 
continuous-time  process  of  this  model.  Therefore,  to  give  more  flexibility  of 
modeling  while  keeping  the  same  properties  as  the  fractional  Brownian  model. 
Fractional  Differencing  model  will  be  introduced  in  the  next  section. 

3.3.  Fractional  Differencing  Model 

3.3.1. 1-D  Fractional  Differencing  Model 

The  1-D  Fractional  Differencing  model  was  suggested  by  Hosking  f  101], 
generalizing  the  well-known  AR1MA  model  of  Box  &  Jenkins  [23],  which  was 
originally  designed  to  model  a  non -stationary  random  process.  A  typical  1-D  first- 
order  Fractional  Differencing  model  is  as  follows. 

_  C 

y(/)  =  (l-z-,)"yC(/)  (3.3.1. 1) 

where  z  is  the  delay  operator  in  the  direction  of  /,  and  £(/)  is  white  Gaussian 
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noise. 

This  random  process  model  has  the  property  of  being  flexible  to  explain  both 
long-term  and  short-term  correlation  structures  of  a  time  series  depending  on  the 
values  of  the  fractional  differencing  parameter  c,  and  it  shares  the  basic  properties 
with  Fractional  Brownian  motion  defined  by  Mandelbrot  [151]. 

Theorem  3.3. 1.1:  The  spectral  density  of  {  y(l)  }  is  s(as)  =  (2sin-i-0))"c  for 
0  <  to  <  k  and  s  (to)  =  as~c  as  co  — »  0,  assuming  that  o^2  =  1. 

_  C 

(Proof)  Let  B  (z)  =  (l-z-1)  2  .  Then,  5(2)  can  be  written  by  B  (z)B~l(z)o^2.  Thus, 
since  o^2  =  1,  we  have  s  (as)  =  B  (eJ(a)B  (e  The  result  follows  on  substitution  of 

_ 

B(eJ(»)  =  (  1 2. 


Theorem  3. 3. 1.2:  The  autocovariance  function  of  (y(l)}  is 


P*  =  2  2  sin(~) - - — -r(l-c),  for  -1  <  ^  <  -±-,  c*0.  (3.3.1.2) 

2  r(*+i~|)  2  2 

and  the  autocorrelation  function  is  given  by 

F(l-y)  nk+~) 

Pi  = - • - 

Hy)  T(k+\-±) 


(3.3.1.3) 
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P*  « 


ro-f) 

r<f) 


(3.3.1.4) 


(Proof)  The  autocovariance,  if  it  exists,  will  be  given  by 


2n 


Costco  s  (co)  d(& 


(3.3.1. 5) 


Jt  c 


=  f  2  2  cos/cto  (sinco/2)  c  dco. 


(3.3.1. 6) 


Using  the  standard  formula, 
ji 

J  sinv-Ix  cosax  dx  = 


7t  cosajt/2 


>v-l 


_ . v+a+1  v-a+1  .  ’ 

v  -fi( — , — ) 


(3.3. 1.7) 


we  can  get  (3. 3. 1.2).  Here,  fl(v)  is  Beta  function.  Also,  from  the  definition  of  the 
autocorrelation  function,  p*  =  — .  Therefore,  (3. 3. 1.3)  follows  on  substitution  of 

and  ft{).  As  &—»«>,  T(k+a)/r(k+b)  can  be  approximated  by  ka~b,  using  the  standard 
approximation  derived  from  Sheppard’s  formula.  Therefore,  it  follows  that  the 
autocorrelation  is  given  by  (3.3. 1.4). 


Definition  3.3.1:  When  c  <  0,  the  Fractional  Differencing  process  (3.3.1 .1)  has  a 
short-term  memory,  and  when  c  >  0,  it  has  a  long-term  memory. 
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From  Theorem  3.3. 1.1  and  3. 3. 1.2,  when  c  >  0,  the  autocorrelation  function  p*  is 
positive  and  decays  monotonically  and  hyperbolically  to  zero  as  the  lag  increases,  and 
the  corresponding  spectral  density  s  (to)  is  concentrated  at  low  frequencies:  s  (co)  is  a 
decreasing  function  of  to  and  s  (co)  -» as  to  -»  0.  Similarly,  when  c  <  0,  s  (to)  is 
dominated  at  high  frequencies:  s(to)  is  a  increasing  function  of  to  and  vanishes  at 
co  =  0.  Therefore,  when  c  <  0  or  c  >  0,  the  Fractional  Differencing  process  (3.3. 1.1) 
has  a  short-term  or  long-term  memory,  respectively. 

Corollary  3.3.1:  The  Fractional  Differencing  process  is  the  discrete  version  of 
fractional  white  noise  process,  and  shares  the  basic  properties  with  Fractional 
Brownian  motion. 


(, Proof)  Brownian  motion  is  a  continuous  time  stochastic  process  B(t)  with 
independent  Gaussian  increments.  Its  derivative  is  the  continuous-time  white  noise 
process,  which  has  constant  spectral  density.  Fractional  Brownian  motion,  BF(t),  is  a 
generalization  of  these  processes.  Then,  Fractional  Brownian  motion  with  parameter 

F ,  usually  0  <  F  <  1,  is  equal  to  the  (y  -  F)th  fractional  derivative  of  Brownian 


motion  in  Riemann-Liouville  sense.  The  continuous-time  fractional  noise  process  is 
then  the  derivative  of  Fractional  Brownian  motion,  thus  it  may  also  be  thought  of  as 


the  (y  -  F)th  fractional  derivative  of  continuous- time  white  noise.  Therefore,  the 


Fractional  Differencing  model  (3. 3. 1.1)  is  the  discrete  version  of  this  continuous-time 
fractional  white  noise  process,  and  it  shares  some  properties  with  Fractional  Brownian 
motion  [101]. 
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3.3.2.  2-D  Fractional  Differencing  Model 

In  3-D  textural  surface  image,  the  fractional  differencing  parameter,  which  is 
‘fractal  scaling  parameter’  in  the  terminology  of  Pentland  [173, 176],  indicates  the 
roughness  of  the  surface,  that  is,  as  the  value  of  the  fractal  scaling  parameter 
increases,  the  model  represents  a  3-D  surface  textured  more  roughly.  However, 
Pentland ’s  model  based  on  2-D  fractional  Brownian  motion  process  with  one  fractal 
scaling  parameter  has  limited  modeling  ability,  because  the  fractal  of  surface  is 
assumed  to  be  spatially  isotropic.  Thus,  Pentland’s  model  is  not  flexible  enough  to 
represent  the  wide  range  of  different  texture  patterns  encountered  in  practice, 
especially  for  the  non-stationary  random  texture.  In  other  words,  this  model  can  tell 
how  rough  the  surface  is,  but  can  not  tell  which  pattern  the  surface  is  covered  with. 
This  limitation  can  be  overcome  by  using  the  2-D  Fractional  Differencing  models  as 
follows. 

3.3.2.I.  First-order  Fractional  Differencing  Model 

As  mentioned  in  previous  sections,  Pentland’s  fractional  Brownian  motion  model 
with  one  fractal  scaling  parameter  is  not  flexible  enough  to  present  the  different 
texture  patterns.  However,  even  from  the  first-order  Fractional  Differencing  model, 
we  could  get  more  flexibility  of  modeling,  because  it  contains  two  different  fractal 
scaling  parameters,  c,  d,  in  directions  of  l\  and  1 2-  This  model  does  not  require  the 
stationarity  of  random  texture.  The  typical  first-order  Fractional  Differencing  model  is 


as  follows. 
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>>(/i,/2)  =  (l-zr1)  2(1  -z2  ')  2^l\Ji)  (3. 3. 2.1.1) 

for/1,/2  =0,1, . N-l  . 

The  corresponding  DFT  of  this  function  is 

-f-  -2  2n^  ~ 

Y(k1,k2)  =  (l~e  N  )  2  (\-e  N  )  2  W(kuk2),  (3.3.2.1.2) 

where  Zj  is  the  delay  operator  associated  with  /,-,  t,(l\,l2)  is  an  i.i.d.  Gaussian 
sequence,  and  W(k\,k2 )  is  the  corresponding  DFT. 

3.3.2.2.  Second-order  Fractional  Differencing  Periodic  Model 


The  typical  second-order  Fractional  Differencing  Periodic  (FDP)  model  is  as 
follows. 

_  C 

y(/1,/2)  =  (l-2^a),zr1+zr2)  2 

_d_ 

*  (l-2c<7ito2  Z2_1+z2  2)  2  £(^1^2)  (3. 3. 2. 2.1) 

for  /1  ,/2  =  0,1,...., N-l  . 

The  corresponding  DFT  of  this  function  is 

..  *1  C 

-) 2k—  -./4k— 

Y(k\,k2)  =  (\-2cosas\e  N+e 

'2k  4k  ^ 

•  (l-2cwco2e~;  ^ )2  W(kuk2),  (33.2.2.2) 

where  z,  is  the  delay  operator  associated  with  /,,  C,(l\,l2)  is  an  i.i.d.  Gaussian 
sequence,  and  W(k\,k2)  is  the  corresponding  DFT. 
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This  model  has  four  different  parameters,  c  and  d  for  the  fractal  scales  and  0)] 
and  o>2  for  the  frequencies  of  pattern  in  the  direction  of  l\  and  1 2,  respectively.  Thus, 
this  model  represents  the  roughness  of  the  surface  and  the  pattern  of  the  texture  image 
at  the  same  time  even  with  the  different  values  for  the  direction  of  l 1 2  separately. 

3.3.3.  Parameter  Estimation 

Because  of  the  flexibility  and  the  simplicity  of  Fractional  Differencing  model,  it 
is  attractive  in  modeling  various  kinds  of  time  series,  including  2-D  random  texture 
image.  However,  the  estimation  of  the  fractal  scaling  parameters  c,  d  and  the 
frequency  parameters  coj,  0)2  is  difficult  and  this  difficulty  further  delayed  the 
application  of  this  model.  There  have  been  several  approaches  for  estimating 
parameters  in  various  kinds  of  fractional  differencing  time  series  since  Hosking 
introduced  Fractional  Differencing  model  [101].  For  example.  Granger  and  Joyeux 
[89]  approximated  this  model  by  a  high-order  auto-regressive  process  and  estimated 
the  differencing  parameter  by  comparing  variances  for  each  different  choice  of  it. 
Lapsa  [121]  suggested  a  maximum-likelihood  estimator  in  the  frequency  domain  and 
showed  the  consistency  of  the  estimator.  This  frequency  domain  analysis  was  further 
studied  by  Eom,  and  a  hybrid  method  of  least-squares  and  maximum-likelihood 
estimations  was  recently  proposed  to  estimate  the  fractal  scaling  parameters  and  the 
frequency  parameters,  respectively  [67].  In  this  thesis,  a  least-squares  or  a 
maximum-likelihood  estimation  algorithm  will  be  applied  to  estimate  the  parameters, 
based  on  Eom’s  algorithm. 
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3.3.3.I.  First-order  Fractional  Differencing  Model 


The  estimation  of  fractal  scaling  parameters  c  and  d  in  the  first-order  fractional 
differencing  model  (3. 3.2.1. 1)  can  be  done  by  a  simple  least-square  estimation  scheme 
in  the  frequency  domain  based  on  a  representation  of  the  logarithm  of  the  process 
which  is  linear  in  the  parameters  as  follows.  By  applying  the  logarithm  operator  to 
(3. 3.2. 1.2),  we  can  obtain 

log \Y(kl  ,k2)  |  =  -ylog  |  l-e~J  n~"  I  -ylog  |  \-e~l2Uir  |  +  log  |  W(k  j  ,k2)  I 


=  --|log|2sin(^-)|  -|log|2sin(-^)|  +log|iy(^1,A:2)| 


(3.3.3.1.1) 


c  Jt^i  d  %k2 

|-log|2sin(— )|  --|log|2sin(-^i)|  -a  +  V(kuk2) 


,  for k\,k2  =  0, 


(3.3.3. 1 .2) 


where  a  =  -£[log|W(£i,/:2)|]  and  V'(&i,£2)  =  log|VF(;t1,Jt2)|  +«. 


Then  9  =  ( c,d ,  a)  can  be  estimated  by  minimizing  the  following  cost  function. 

N_  i  N/2  ^ 

J (0,o)1,co2)  =  £  £  (log |  Y {k\,k2 1  +  Iog|2sin-^— 1 

*,= o  Jfc2= o  1  N 


d  kk2  ~ 

+  —  log|2sin-^-|  +  a)2 


(3.3.3. 1.3) 
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N- 1  M 

=  I  I  (log|y(*,.*2)| -0rs«l.*2))2  (3.3.3.1.4) 

*,=0  *2=0 


Here, 


Q(ki,k2)  = 


-1  n  kx 

— -  log  1 2sin— — —  | 
2  N 

-1  7t£2 

—  log  1 2sin— — —  | 
2  A' 

-1 


(3.3.3.1.5) 


Thus,  the  estimated  values  will  be 


- . .  r  *-i  M  T  ,  *-i  H 

[c,d,af  =  (Z  £  Q(ki,k2)QT(ki,k2))~H  £  £  g(*1,*2)log|y(*1,*2)l) 

*1=0  *2=0  *,=0  *2=0 


(3.3.3.1.6) 


3.3.3.2.  Second-order  Fractional  Differencing  Periodic  Model 

The  estimation  of  parameters  in  the  Fractional  Differencing  Periodic  model 
(3. 3. 2. 2.1)  can  be  done  in  the  frequency  domain  [127]. 

For  estimation,  all  parameters  can  be  estimated  directly  from  the  given  data 
Y(k\,k2),  if  we  can  obtain  the  likelihood  function  of  \Y{k\,k2)\.  Then,  since  the 
noise  sequence  i \{l\,l2 )  is  assumed  to  be  white  Gaussian,  W(k\,k2)  and  Y(k\,k 2) 
follow  the  Rayleigh  distribution  as  in  the  following  theorem. 


Theorem  3.3.3. 1 :  The  modulus  of  the  DFT  of  the  noise  sequence,  {\W(k\,k2)  \ , 


ky=0,l  ,  k2=0,l 


N/2 


}  and  {\Y(kuk2)\,  k  i=0,I  , 
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k2=0,l,..., 


N/2 


}  are  the  white  sequences  with  the  following  Rayleigh  densities. 


f \W(kuk2)\(W)  = 


2 W  ,  -W2 
exp  { 


pN 2  p N‘ 

0 


,  W>  0 
,  otherwise 


f \Y(ki.k2)\  00  ~  1 


2s\,t2nkuk2)  -s2tltkj2(k  1,*2) 

—  exp/ - '-2-^ - } 


pNd 


pN‘ 


0 


,  y(*i,*2)*o 

,  otherwise 


where 


2nk\  2  2kIc2 

skuk2  =  1 2cos( — — — )  -  2coscox  |  1 2cos( — — — )  -  2cosa>2 1 


(. Proof)  Let  W(kx,k2)  be  the  DFT  of  a  white  noise  sequence  £(/ \,l2)  with  a  normal 
distribution,  that  is, 

C(/,,/2)~mp)  ,for  l ul 2  =  1,2,. ~,N  (3.3.3.2.1) 

Then,  W(kx,k2)  will  follow  another  normal  distribution  as  follows  [27J. 

W(k x  ,k2)  -  N(N'E  &(lul2)],pN2)  (33.3.2.2) 

~N(0,pN2)  (3.3.3.23) 

and 

Rc(W(kx,k2)),  lm(W(kx,k2))~N(0,^-)  (333.2.4) 

where  Re(W)  and  Im(VV')  are  the  real  and  the  imaginary  parts  of  W(k  x,k2), 
respectively. 
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Now,  consider  the  density  function  of  \W(kx,k2)\.  Since  W(kx,k2)  is  a 
complex,  defining  W  be  W(k\,k2), 


\W(kuk2)\  =  ^Rc2(W)  +  lm2(W) 


(33.3.2.5) 


Define  0  =  tan  1  (7^77—).  Then,  we  can  have  the  following  relations. 
Re(W0 


Re(WO=  \  W{k\,k2)  \  cos0  (33.3.2.6) 

Im(W)  =  |  W  (k\,k2)  |  sin0  ,for\W\  >0  (333.2.7) 

Thus,  the  joint  probability  density  of  |  W  |  and  0  can  be  obtained  from  the  Jacobian 
and  the  joint  probability  density  of  Re(W0  and  Im(W)  as  follows. 

/|W|  <,(1  W|,<»  =  fRe(W)lm(W)(Re(W)M(W))  ,for\W\  >0  (333.2.8) 

=  l^l/Re(W)(|W/ko.v<S))/Im(VV)(|W'|i-m<J))  (33.3.2.9) 


-(Re2(Wy+lm2(W)) 

_  JiLLe  p^2 

TtpN2 


(3.3.3.2.10) 


where  Jacobian  = 


ai  w  l 
dRe  ( W ) 
90 

dRe  (W) 


a  l  w  1 

d!m(W) 

90 

dim  (W') 


(3.3.3.2.11) 


1 

\W\ 


(3.3.3.2.12) 


Iherefore,  since  |W|  and  0  are  independent,  and  the  density  function  of  0, /^(0),  is 
1 

2tc  ’ 


f\  W I ,  *(  I  w  1 , 9)  =  / 1  w  1  ( |  w  I  )-/0(0) 


(3.3.3.2.13) 


and 
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2W  - W 2 

f\w\m  =  J,pN2ex^  pN2>  W>Q  (3.3.3.2.13) 

0  ,  otherwise 

Now,  set 

c  d 

2nk\  T  2nk2  T 

Sk  k ,  =  \2cos( - )-2cas'C01|  |2 cos( — - — )  -  2cosg>2  I  (3.3.3.2.14) 

N  N 

Then,  we  can  have 

|r(*i,*2)l=*Jt,.*r1  \w&iM)\  (3.3.3.2.15) 

Thus,  the  density  function  of  |  K (Jfe  i , A: 2 )  I  can  be  represented  by  following 
equations. 

f  \Y(kltk2)\  (Y)  =  skltk2  f  |W(*,.*2)|  (Skuk2Y(k  1  ,£2))  (3.3.3.2.16) 

2s2t„t,Y(kl,k2)  , 

pN2  eXPf  pN2  ;  ,  K(*!,*2)>0 

3  ,  otherwise 

U 

From  these  probabilistic  properties,  the  estimation  of  parameters  can  be  done  by 
a  hybrid  method  of  Least  Square  and  Maximum  Likelihood  estimations,  which  was 
suggested  by  Eom  [67],  For  LS  estimation,  if  the  values  of  parameters  a>i ,  0)2  are  set, 
then  0  =  (c\d,  a)1  can  be  estimated  by  minimizing  the  following  cost  function. 

n-i  h2J  c  2  Kkx 

J  (6,C0i  ,to2)  =  X  £  (log  |  K  (/c  1,^2 1  +  —  log  1  — — - 2cos®\  | 

*,=o  *2=o  1  N 
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d  , 

+  —  log  1 2 cos—jj - 2c os (£>2 1  +  a) 


(3.3.3.2.17) 


N-l  H 

:  £  £  (logiy^!,^)! -QTQ(kl,k2))2 

*,=  0  *2=0 


(3.3.3.2.18) 


where  a  =  -£[log|  W(A: i ,^2) I ]»  0  =  (c,d, a)7. 


Here, 


-1 

—  log  1 2 cos—jj - 2cos(Oi  | 

Q(k  i  ,k2)  =  -l  2tc*2  .  . 

—  log  1 2 cos— - 2cos(02  | 


(3.3.3.2.19) 


Thus,  the  estimated  values  will  be 

yv-i  [Nl2\  N-i  M 

[c,J,cx]7'  =  (  X  £  Q(kl7k2)QT(/ci,k2)r\Z  I  Q0cuk2)log\Y(k,,k2)\ 
*1=0  *2=0  *,=0  *2=0 


(3.3.3.2.20) 


Also,  ML  estimators  of  tOj ,  0)2  can  be  calculated  by  maximizing  the  log- 


likelihood  function  L(Y ;0,(Di  ,o>2)  with  a  estimated  from  above. 


/-!  K 


h’-l  [Nl\ 


N-i  l  j  N-i  r*\  on2 

L(YAc b,,<o2)=  £  £  log|r(*1,*2)|-  £  £  log(-^— ) 

*,=0  *2=0  *,=0  *2=0  l 


c  N-\  2nk{  M  2nk2 

+  —N  £  log  1 2 (cos—— - coscoj )  1  +dN  £  log  1 2{cos— - coso>2)  | 

1  *,=0  "  *2=0  " 
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1  AM  ru\  2%k\  c  2nk2  d  , 

- £  £  | 2(c<?5 — — - COSO)! ) I  •  1 2 (cos—— - coso^)!  *  |  Y(k {,k2)  \ 2 

P N2  k%  N  N 

(3.3.3.2.21) 

where  p  is  a  variance  of  £(/ \,l2)  and  can  be  estimated  by  the  following  equation 
in  mean  square  sense. 

p  =  — ^-exp  fy- 2a — ~t}  (3.3.3.2.22) 

N2  6N2 

where  y  is  Euler’s  constant  (=  0.5772157)  [67]. 

Therefore,  the  estimation  scheme  can  be  summarized  as  follows: 

(Estimation  Algorithm) 

Step  1 :  Choose  resonant  frequencies  0)! ,  o>2  in  the  range  of  [0. yl- 

Step  2:  With  the  given  values  of  Wj  and  o>2,  estimate  c,  d,  and  a  by  LS  estimation 
algorithm  (3.3.3.2.20). 

Step  3:  Using  a,  compute  the  estimate  of  the  variance  of  p,  £,(li,l2),  by  equation 
(3.3.3.2.22). 

Step  4:  Using  the  estimates  c,  d  and  p  found  in  Steps  2  and  4,  maximize  the 
likelihood  function  given  by  (3.3.3.2.21)  with  respect  to  coj  and  o>2- 

A  A 

Step  5:  Using  the  estimates  CDi  and  co2,  repeat  Step  2  to  Step  4  until  the  estimates 
have  no  significant  change  in  successive  iterations. 
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3.4.  Conclusions 

In  this  chapter,  various  types  of  Fractional  Differencing  model  have  been 
considered  and  compared  with  the  Fractional  Brownian  model.  The  Fractional 
Brownian  model  has  been  used  for  measuring  the  roughness  of  a  surface  by 
researchers.  However,  the  modeling  ability  of  Fractional  Brownian  motion  noise  is 
limited,  because  this  model  has  only  three  variable  parameters,  mean,  variance  and 
fractal  scale  F,  thus,  those  are  not  enough  to  model  the  wide  range  of  randomness 
encounted  in  practice.  Also,  it  is  relatively  difficult  to  estimate  the  fractal  scale  F  and 
to  synthesize  with  it,  because  of  its  continuous-time  process.  To  give  more  flexibility 
of  modeling  keeping  the  same  properties  of  Fractional  Brownian  motion  process,  the 
Fractional  Differencing  model  is  suggested.  The  Fractional  Differencing  model  is  a 
discrete  version  of  the  Fractional  Brownian  model.  Thus,  this  model  shares  the  same 
properties  of  the  Fractional  Brownian  process  and  has  a  relatively  simple  structure.  In 
this  model.  Fractal  scale  can  have  two  different  values  in  the  directions  of  x  and  y,  and 
this  gives  more  flexibility  for  modeling  non-isotropic  distributed  random  texture. 
Additionally,  this  model  can  be  extended  to  the  second-order  of  process  to  provide 
another  pair  of  frequency  parameters  CD],  o>2.  Using  this  2nd-order  Fractional 
Differencing  model,  the  roughness  and  the  texture  pattern  of  a  surface  can  be 
represented  simultaneously.  Parameter  estimation  schemes  for  each  model  are  also 
developed  based  on  the  least-square  estimation  and  the  maximum  likelihood 
estimation. 
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CHAPTER  4 

SHAPE  FROM  A  SHADED  AND  TEXTURED 
SURFACE  IMAGE 


4.1.  Introduction 

An  important  task  in  computer  vision  is  the  recovery  of  3-D  scene  information 
from  single  2-D  images.  3-D  analysis  of  an  image  can  be  broken  down  into  two  main 
categories,  Shape-from-shading  and  Shape-from-texture.  The  Shape-from-shading 
technique  uses  the  reflectance  map  which  shows  scene  radiance  as  a  function  of  the 
surface  gradient  and  the  distribution  of  light  sources  to  extract  3-D  surface 
information  from  image  data  [99, 144, 172],  On  the  other  hand,  the  Shape-from- 
texture  analysis  technique  uses  the  texture  pattern  instead  of  shading  to  extract  3-D 
structure.  Since  texture  gradients  behave  like  intensity  gradients,  the  shape  of  a 
surface  can  be  inferred  from  the  pattern  of  a  texture  on  the  surface  by  applying 
statistical  texture  analysis  [6, 135, 206, 221  \. 

However,  for  describing  a  natural  scene  image,  each  of  the  above  approaches  has 
its  own  limitation.  The  Shape-from-shading  technique  is  applicable  only  under  the 
assumption  that  the  surface  is  smooth  enough  to  have  clear  radiance  information, 
while  the  Shape-from-texture  technique  requires  the  surface  to  be  relatively  complex 
so  that  texture  information  can  be  extracted.  Thus,  neither  technique  is  suitable  to 
recover  3-D  structure  information  from  a  natural  scene,  because  both  radiance  and 
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texture  information  coexist  within  the  surface  of  a  natural  scene.  Therefore,  a  robust 
technique  is  needed  to  handle  this  shortcoming.  Recently,  the  fractal  scaling 
parameter  was  introduced  to  measure  the  coarseness  of  the  surface,  and  applied  to 
represent  the  natural  scene  surface  [176],  However,  this  fractal  model  is  not  enough 
to  represent  the  real  3-D  texture  image,  because  even  though  two  surfaces  are 
estimated  to  have  the  same  fractal  scales,  these  surfaces  can  have  different  texture 
patterns. 

In  this  chapter,  a  composite  model  of  Shape-from-shading  and  Shape-from- 
texture  is  developed  to  represent  a  3-D  surface  image  considering  the  scene  image  as 
the  superposition  of  a  smooth  shaded  image  and  a  random  texture  image.  The 
orthographical  projection  is  adapted  to  take  care  of  the  non-isotropic  distribution 
function  due  to  the  slant  and  tilt  of  a  3-D  texture  surface,  and  the  Fractional 
Differencing  Periodic  model  is  chosen  because  this  model  is  able  to  simultaneously 
represent  the  coarseness  and  the  pattern  of  the  3-D  texture  surface  with  the  fractional 
differencing  parameters  c,  d  and  the  frequency  parameters  coj,  (1)2.  and  it  has  the 
property  of  being  flexible  enough  to  synthesize  both  long-term  and  short-term 
correlation  structures  of  random  texture  depending  on  the  values  of  the  fractional 
differencing  parameter  c  and  d.  Since  the  object  is  described  by  a  model  involving 
several  free  parameters  and  the  values  of  these  parameters  are  determined  directly 
from  its  projected  image,  it  is  possible  to  extract  3-D  information  and  texture  pattern 
directly  from  the  given  intensity  values  of  the  image  without  any  pre-processing. 
Thus,  the  cumulative  error  obtained  from  each  pre-processing  can  be  minimized.  For 
estimating  the  parameters,  a  hybrid  method  which  uses  both  the  least  square  and  the 
maximum  likelihood  estimates  is  applied  and  the  estimation  and  the  synthesis  are 
done  in  frequency  domain  based  on  the  local  patch  analysis.  By  using  this  model,  the 
integrability  problem  which  might  occur  in  spatial  domain  analysis  can  be  avoided, 
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because  only  one  inverse  Fourier  transform  needs  to  be  taken  at  the  end  of  procedure 
to  get  the  whole  image. 

The  organization  of  this  chapter  is  as  follows.  In  section  4.2.,  we  introduce  the 
image  model  i(l\,l2)  which  is  obtained  by  superposing  the  deterministic  function 
x(l 1 2)  and  the  random  function  y(l\,l2)>  and  the  relationship  between  the  different 
directions  of  a  3-D  surface.  In  section  4.2. 1-3,  an  estimation  scheme  for  the 
illumination  direction,  the  modified  reflectance  map  function  for  cc  (/ 1 , / 2 ) »  and  the 
orthographically  projected  Fractional  Differencing  Periodic  function  for y{l\Ji)  are 
introduced.  Section  4.3.1  outlines  the  estimation  scheme  for  the  parameters  in  the 
composite  model.  Section  4.3.2  then  discusses  some  simulation  results  carried  out  to 
demonstrate  the  performance  of  the  proposed  algorithm,  followed  by  Section  4.4 
which  concludes  the  paper. 

4.2.  Model  Of  3-D  Texture  Surface  Image 

The  surface  shape  is  usually  defined  in  terms  of  the  viewer’s  coordinate  system. 
This  system  has  axes  / 1 ,  lj,  1 3  with  the  / 3  axis  in  the  viewing  direction.  The  observed 
intensity  function  <(/ 1,/2)  of  local  shape  of  a  3-D  surface  can  be  considered  as  the 
sum  of  a  deterministic  function  and  a  statistical  random  function  y(l\,l 2), 

whose  expected  value  is  zero: 

i(l  1  Ji)  -  x0 1 J2)  +  y(l  1  X2)  (4.2.1) 

and 

E[I(/,,/2)1=a(/1,/2)  (4.2.2) 

Thus,  x(l  1 ,/ 2)  can  be  simply  estimated  by  smoothing,  i.e.,  taking  the  average  of 
intensity  values  in  the  proper  size  of  window.  Then  y{l j,/2)  can  be  estimated  by 
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subtracting  this  function  from  the  original  image  / i , / 2 )-  Figure  4.1  shows  this 
superposition  in  a  simple  1-D  case.  Note  that  these  intensity  functions  do  not 
represent  the  actual  shape  of  image,  because  the  shading  variations  are  caused  by 
changes  in  surface  orientation  relative  to  the  illumination  direction  L.  Therefore, 
x(l  1 , / 2)  and  y(l\,l2)  need  to  be  projected  to  the  illumination  direction  to  extract  the 
shape  information 

Let  <3i,  be  the  slant  and  tilt  of  illumination  direction  L  in  the  coordination 
system  with  l\,  1 2,  h.  The  direction  L  induces  a  coordinate  system  with  axes 
i  [  t'2,  i 3 ,  and  this  is  derived  from  the  following  transformation. 

i  1  coso^cosx^  cosa^sint/.  -sinx^  l\ 

i  2  =  ~sinxL  cosX£.  0  l2  (4.2.3) 

i  3  sino^cosx^  sino^sinx^  cosg^  /  3 

Here,  /  3  is  the  illumination  direction  L,  and  /  3  is  the  viewing  direction. 

The  surface  normal,  N,  is  another  important  direction  of  3-D  geometry.  If  we 
define  a  new  coordinate  system  with  axes,  mx  m2,  m3,  having  the  ^3-axis  in  the 
surface  normal  direction,  these  axes  can  be  derived  from  the  coordinate  system  of 
viewing  direction,  i.e.,  l\,  1 2,  / 3,  by  using  a  coordination  transform  similar  to  the 
above  with  different  values  of  tilt  x  and  slant  a. 

^  1  cosacosx  cosasinx  -sinx  m\ 

1 2  =  -sinx  cosx  0  m2  (4.2.4) 

[3  sinocosx  sinasinx  coso  m3 

Refer  Figure  2.1  for  the  relations  between  these  three  different  directions  and  the  four 
different  tilt  and  slant  p;irameters,  a 1,  x^,  G,  and  x. 

Let  x\m\,ni2)  and  >’'(mi-ffl2)  be  a  deterministic  and  a  random  function 
respectively,  defined  on  the  surface  normal  image  plane,  i.e.,  m 3=0.  Then, 


deterministic  functions  in  1-D  case. 


z'(ml,m2),  defined  as  the  sum  of  these  two  functions,  is  merely  the  depth  function  z  , 
observed  from  the  surface  normal.  However,  note  that  y'(m1,m2)  is  not  the  rotated 
function  of  >’(/,, /2)  but  the  projected  function  ofy(/,,/2)  to  the  m,-m2  plane  so  that 
they  may  satisfy  the  superposition  property.  Figure  4.2  depicts  these  relations 
between  functions. 

4.2.1.  Model  Of  The  Deterministic  Component  jt 


As  discussed  before,  a  3-D  surface  image  can  be  considered  to  be  the 
superposition  of  a  texture  image  and  a  smooth  shaded  image.  This  smooth  shaded 
image  can  be  represented  by  the  deterministic  function  j t(/1,/2).  Thus,  if  the 
illumination  direction  values  for  the  whole  image  are  given,  the  surface  orientation 
parameters  slant,  o,  and  tilt,  x,  can  be  estimated  from  x{l  \,l2)  representing  a  smooth 
surface  by  Shape-from-shading  analysis. 

Pioneering  work  on  the  inference  of  shape-from-shading  was  done  by  Horn  [99] 
and  his  co-workers.  To  extract  the  3-D  shape  function  //(•,•)  from  a  single  2-D  image, 
they  used  the  reflectance  map,  which  shows  the  intensity  of  the  image  as  a  function  of 
the  surface  gradient  and  the  illumination  direction. 


x{l\J2) 


pcosT^sino^+r/sint^sina^+coso^ 

(p2+q2+])V2 


(4.2.2. 1) 


where 


/?(p,q)  :  Reflectance  map  function 

p  =  4;,ni'j2)’  q  =  jqHU'’h) 


H(l  i  ,/2) :  3-D  shape  function  from  the  viewing  direction. 
:  Tilt,  slant  of  the  illumination  direction 
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Figure  4.2:  3-D  Geometry  of  the  functions  x,  y ,  x\ y'. 


Here,  from  the  relationship  between  the  tilt  x,  the  slant  o  of  the  surface  and  p,  q 


X  =  tan  x{— ), 
P 


a  =  cos  —  — 

']p2+q2+ 1 


p  =  tanocosx,  q  =  tancsinx  , 

we  can  modify  (4.2.2. 1)  to  the  function  of  o,  x,  ct^,  and  X£,  as 

sinacosxcosX£,sin(j£,  +  sinosinxsinx^sinc/,  +  cosacosa/. 


(4.2.2. 2.a) 

(4.2.2.2.b) 

(4.2.2.3) 


Construction  of  3-D  shape  can  be  achieved  by  solving  a,  x  in  terms  of  x(li  J2)  at 
each  point,  and  integrating  those  values.  However,  this  approach  needs  the  solutions 
of  at  least  two  difficult  problems.  First,  this  is  an  ill-posed  problem  because  there  are 
two  unknown  parameters  and  only  one  equation  to  be  solved.  Thus,  we  need  an 
additional  constraint  to  have  a  unique  solution.  Second,  the  final  integrated  shape  can 
be  different  from  the  original  shape,  due  to  the  cumulation  of  estimation  errors. 

To  get  a  unique  solution,  the  calculus  of  variation  methods  were  used  by 
minimizing  the  estimation  error  after  adding  one  constraint  for  the  smoothness  [100] 
or  integrability  [76].  To  handle  the  cumulative  error,  Pentland  [172]  suggested  local 
shape  analysis  which  deals  with  only  the  local  areas  instead  of  a  whole  image. 
However,  Pentland’s  technique  has  severe  trouble  in  integrating  all  local  area  surface 
information.  Recently,  Pentland  [169]  developed  another  technique  for  solving  the 
integrability  problem.  He  suggested  analysis  in  the  frequency  domain,  instead  of  in 
the  spatial  domain.  By  using  this  method,  the  integrability  problem  can  be  avoided, 
because  only  one  inverse  Fourier  transform  needs  to  be  taken  at  the  end  of  the 
procedure.  However,  since  the  calculation  of  a  convolution  is  required  in  the 
frequency  domain  to  handle  a  simple  multiplication  operation  in  the  spatial  domain, 
calculation  will  be  complicated.  In  this  paper,  these  ill-posed  and  the  integrability 
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problems  will  be  handled  by  an  additional  constraint  from  the  texture  pattern  and 
frequency  domain  analysis,  respectively.  (This  will  be  discussed  in  detail  later.) 

4.2.2.  Model  Of  The  Random  Component  y 

The  random  component  3^ (/ 1 » / 2 )  of  the  intensity  function  /(/ 1 ,/ 2)  can  not  he 
simply  obtained  by  taking  coordinate  transformation  to  2),  defined  on  the 

surface  normal  plane  image.  Because  the  random  function  y(l  j  ,/2)  was  assumed  to  be 
a  2-D  random  function,  whose  expectation  value,  E[y  (/ 1 ,/ 2)J»  >s  zero>  to  satisfy  the 
superposition  properties,  this  function  should  rather  be  considered  as  the  function 
obtained  after  projecting y\tn\,m 2)  to  the  viewer’s  direction  (as  shown  in  Figure  4.2). 
Thus,  if  we  can  model  y'{m\,mi)  properly,  >’(/], I2)  will  be  obtained  by  projecting 
this  to  the  viewer’s  direction  L. 

4.2.2. 1.  Fractional  Differencing  Periodic  Model  Fory'(mi,m2) 

The  random  function  y'(m  \  .wj),  defined  on  the  surface  normal  plane  image,  can 
be  approximated  by  a  2-D  random  field  model,  which  is  distributed  over  the  surface 
normal  plane  (as  shown  in  Figure  4.1).  Note  that  since  this  model  is  based  on  a  2-D 
texture  model,  it  will  fit  the  plane  surface  more  than  any  other  shape  of  surface.  Thus, 
we  apply  this  model  to  the  local  shape  of  the  image  because  the  local  shape  will  be 
closer  and  closer  to  the  plane  surface  when  we  take  a  smaller  patch. 

Among  the  various  random  field  models,  the  Fractional  Differencing  Periodic 
model  is  chosen  to  represent  this  random  texture  surface  in  this  chapter.  This  random 
process  model  has  the  property  of  being  flexible  enough  to  explain  both  long-term  and 
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short-term  correlation  structures  of  a  time  series  depending  on  the  values  of  the 
fractional  differencing  parameter  c,  that  is,  when  c  <  0,  the  fractional  differencing 
process  has  a  short-term  memory,  and  when  c  >  0,  it  has  a  long-term  memory 
(Theorem  3.3.1).  Also  it  shares  the  basic  properties  with  Fractional  Brownian  motion 
defined  by  Mandelbrot  [151]  (Theorem  3.3.2). 

In  3-D  textural  surface  image,  the  fractional  differencing  parameter,  which  is 
‘fractal  scaling  parameter’  in  the  terminology  of  Pentland  [173,176],  indicates  the 
roughness  of  the  surface,  that  is,  as  the  value  of  the  fractal  scaling  parameter 
increases,  the  model  represents  a  3-D  surface  textured  more  roughly.  However, 
Pentland’s  2-D  model  with  one  fractal  scaling  parameter  has  limited  modeling  ability, 
because  the  fractal  of  surface  is  assumed  to  be  spatially  isotropic.  Thus,  Pentland  s 
model  is  not  flexible  enough  to  represent  the  wide  range  of  different  texture  patterns 
encountered  in  practice,  especially  non-stationary  random  texture.  In  other  words,  this 
model  can  tell  how  rough  the  surface  is,  but  can  not  tell  which  pattern  the  surface  is 
covered  with.  This  limitation  can  be  overcome  by  using  the  2-D  Fractional 
Differencing  Periodic  model  as  follows  [67 j. 

C 

y'(m  \  ,nii)  =  (l-2co.vo)| z'r'+  z'\~2)  2 

_d_ 

•  (l-2cavo>i  z’2 ~'+zV2)  (4. 2. 3. 1.1) 
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-jin—  -j4n- f 

'(X-lcosv^e  N  +e  N  )  2  W'(kuk2),  (4.2.3. 1.2) 

where  z',  is  the  delay  operator  associated  with  m,,  ^'(ff!1,/n2)  is  an  i.i.d. 
Gaussian  sequence,  and  \V\k\,k.2)  is  the  corresponding  DFT. 

This  model  has  four  different  parameters,  c  and  d  for  the  fractal  scales  and  0^ 
and  to 2  for  the  frequencies  of  pattern  in  the  direction  of  nt\  and  m2,  respectively. 
Thus,  this  model  represents  the  roughness  of  the  surface  and  the  pattern  of  the  texture 
image  at  the  same  time  even  with  the  different  values  for  the  direction  of  n%\ ,  m2 
separately. 

4.2. 2.2  Orthographically  Projected  Fractional  Differencing  Periodic  Model  For 

> di,i2) 


The  statistical  model  of  the  intensity  function  >’ (/ 1 » / 2 )  cannot  be  simply  obtained 
by  rotating  the  coordinate  axes  because  the  expected  values  of  both  /(mi,ffi2)  and 
v ( / 1 . / 2 )  must  be  zero  over  the  planes  to  satisfy  the  superposition  properties. 
Therefore,  such  a  function  y{l \,lj)  that  satisfies  this  requirement  can  be  obtained  by 
projecting  function  /(m^mi)  orthographically  to  the  viewer’s  image  plane.  Figure 
4.2  shows  this  projection. 

A  new  coordinate  system  of  the  orthographically  projected  image  from  the 
viewing  direction,  /  j  -/ 2 ,  can  be  obtained  from  the  following  two  coordinate 


transformations. 
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>  1' 

cosx  -sinx 

m  1 

/'  2 

sinx  cosx 

m  2 

(4.2.3.2.1) 


and 


r  1 

l\ 

1  0 

w 

h 

0  coso 

I'l 

(4.2. 3. 2. 2) 


Here,  x  is  the  angle  between  m\  and  l'\  axes  and  o  is  the  rotational  angle  based  on 


l'x  -axis  (Refer  Figure  2. 10). 

Hence,  the  coordinate  transformation  of  the  orthographic  projection  between  the 
m  i  -m2  system  and  / 1  -l 2  system  can  be  given  as  follows  [110]. 


h 

h 


cost  -sinx 
sinx  cosx 


■ 

1  0 
0  coso 


r 

-1 

• 

cosx  -sinx 

mi 

sinx 

cosx 

m  2 

L  J 

. 

(4.2. 3.2. 3) 


cos2x+cososin2x  ( l-coso)sinxcosx 
( 1  -cosa) sin xeosx  sin2x+coscrcos2x 


mi 

in  2 


(4.2. 3.2.4) 


Thus, 


*1 

nil 

1 

sin2x+cosocos2x 

(coso-l)sinxcosx 

h' 

mi 

COSO 

(coso-l)sinxcosx 

cos2x+cososin2x 

A  more  detailed  discussion  on  the  orthographical  projection  and  some  examples 
to  demonstrate  this  projection  were  given  in  section  2. 3. 1.1  (Figure  2.1 1). 

As  a  result  of  the  coordinate  transform  (4.2. 3.2. 5),  the  model  of  intensity 
function  y(/ 1,/ 2)  can  be  obtained  from  the  Fractional  Differencing  Periodic  model  of 
y'(m  1  ,m  2)  as  follows. 

Let 
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>’(/i,/2)  =  (l-2costo1z1  '+  z j  2)  2  (1-2cosco2  z2  x+z2  2 )  2£>Q\,h)  (4.2.3. 2.6) 
where  z ! ,  z2  are  the  delay  operators,  corresponding  to  /  j  and  / 2>  respectively. 


By  the  definition  of  DFT,  IV (X:  | 2 )  corresponding  to  white  noise  sequence  C(^  1*^2) 


n-  1  N-i  2n 

W(kuk2)  =  X  Z  ew\-j-jj-(mlkl+m2k2)]  (4.2.3.2.7) 

m,=0/n2=0 

and,  from  (4.2. 3.2.5), 


N- 1  yv-l  2n 

W (k\,k2)=  £  - [((si^t  +  cosocos2!)/:! 

/,=o/2=o  ^coso 


+  (cosa-  l)sin'CcosT&2Vi  +((cosa-  l)sinTcosTfc] 


+  (cos2x  +  cososirrt)&2)(  2)) 


(4.23.2.8) 


Thus,  as  in  1 1 1 1 1,  we  can  define 


sin"’ 1  +  cosocos2!  (coso  -  l)sinxcosx 

,  _  '  cose  '  COSO 

z  1  -  z  1  z  1 


(4.2.3.2.9a) 


(cosO  -  1  JsinXcosx  cos' x  +  cososin2X 


. '  COSO 

Z  -i 


(4.2.3.2.9b) 


n  1  = - |(sitrT  +  cosocos"T)&  1  +  (coso-  l)sinxcostit2],  (4.2.3.2.10a) 

coso 


n2  = - [(coso  -  DsinTcosU]  +  (cos2T  +  cososin2T)&->|  (4.2.3.2.10b) 

coso 


Thus,  from  the  orthographical  projection,  we  can  set 
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W(kl,k2)  =  W'(nl,n2)  (4.2.3.2.11) 

Therefore,  with  these  relations,  we  can  have  the  projected  version  of  Y'(ki,k 2), 
which  is  DFT  of  y(l\,l2 )  as  follows. 

Y(ki,k2)  =  (l-2coso)ie  N  +e  N)2 

-j2n—  -jAn— 

•(l-2cosa>2<T;  71  N  +e  1  K  N  )  2  W(kuk2)  (4.2.3.2.12) 

Note  that  the  fractal  scaling  parameters  c  and  d  remain  same  as  the  ones  before 
projecting  because  of  their  scaling  invariance  property  [101, 125). 

4.3.  Projected  Texture  On  The  3-D  Surface 

As  discussed  in  previous  chapters,  the  intensity  function  of  an  image,  i(l\,l2), 
can  be  represented  by  the  superposition  of  a  deterministic  function  x (l\,l2)  and  a 
random  function  y(l\,l  2).  Therefore,  if  we  can  estimate  all  parameters,  that  is,  a,  x  for 
the  surface  orientation  and  c,  d,  0)j,  0)2  for  the  pattern  of  texture  from  the  intensity 
function  z  directly,  then  we  can  get  better  estimates  than  the  ones  from  the  separate 
procedures  for  each  x  and  >’  function.  It  is  obvious  that  the  estimation  error  from  one 
procedure  will  cause  another  estimation  error,  thus,  the  errors  will  be  cumulated.  For 
that  reason,  in  this  chapter,  a  composite  model  of  Shape-from-shading  and  Shape- 
from-texture  will  be  discussed. 

Consider  the  reflectance  map  function  (4.2.2. 3)  for  x  and  the  2-D 
onhographically  projected  Fractional  Differencing  Periodic  model  (4. 2. 3. 2. 6)  for  y. 
Then  the  intensity  function  z  can  be  represented  by 
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/(/  i  ,/2)  =  etsinGcosxcosx^sinG^  +  sinGsinxsinx^sino/,  +  cosocosg^J 

_ c_  _ d_ 

+  (1— 2cosc0i zj1  +zj2 )  2  (1-2cosco2221+Z22)  2W\,h)  (4.3.1) 

where  8  is  the  normalization  factor,  and  z\  and  z2  are  same  as  defined  by 
(4.2.3.2.9a)  and  (4.2.3.2.9b),  respectively. 

Here,  notice  that  the  normalization  factor  e  should  be  multiplied  to  the 
reflectance  map  function  R(c,x),  because  /?(g,x)  is  the  normalized  function  which  has 
maximum  value  1.  Because  our  model  is  based  on  the  assumption  that  each  patch  is  a 
slanted  and  tilted  plane,  the  value  of  the  deterministic  function  x(l\  ,/2)  is  a  constant 
over  each  patch,  the  best  guess  of  the  value  of  e  is  the  maximum  value  among  the 
average  intensities  of  the  patches, 

i  N'  N’ 

8=  Maximum  (— £  £  /,(/i,/2)]  (4.3.2) 

N  /,=!  ;2=1 

where,  M  :  Total  number  of  patches  in  the  whole  image. 

N' :  Size  of  the  patch 

/,(•,•) :  Intensity  function  of  the  i-th  patch. 

Therefore,  the  value  of  e  can  be  estimated  from  the  whole  image,  before  the  actual 
procedure  on  each  patch.  Thus,  we  assume  e  to  be  given,  just  as  the  illumination 
directions  Ol,Tl. 

The  DFT  of  (4.3.1)  will  be 

I  {k  \  ,l<2)  =  X  (k  \  ,k2)  +  Y(k]  ,k2)  (4.3.3) 

=  /V2e|sinocosTcosT£,sinG£  +  sinGsinxsinx^sinG^  +  cosgcosG/J5(£  j  ,&2) 
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+  (1— 2coscoic 


n  i 

~jln~N 


. ,  n\  c 
-j4n~N  — 


+e 


)  (1-2coso>2^ 


-j  2n 


«2 

77" 


.  . 

-;47t^r  -o- 


+e 


A/ 


)  *W(kuk2) 


(4.3.4) 


where  ti\,n2  are  same  as  defined  by  (4.2.3.2.10a),  (4.2.3.2.10b),  and 


5(*i,*2) 


|l  k  ]  ,k2=  0 
|o  otherwise 


(4.3.5) 


4.3.1  Estimation  Of  The  Parameters  c,  d,  CD] ,  co2,  a,  t 

The  estimation  of  parameters  in  this  projected  Fractional  Differencing  Periodic 
model  (4.3.4)  can  be  done  in  the  frequency  domain,  modifying  the  techniques 
suggested  in  [  127].  This  frequency  domain  analysis  has  an  advantage  over  the  spatial 
domain  analysis.  Since  we  need  only  to  take  the  inverse  Fourier  transform  to  get  the 
whole  image  at  the  end  of  procedures,  we  can  avoid  the  integrability  problem  which 
might  occur  in  spatial  domain  analysis  when  the  estimated  surface  orientation  from 
each  local  shape  has  an  error. 

For  estimation,  all  parameters  can  be  estimated  directly  from  the  given  data 
l{k\,k2),  if  we  can  obtain  the  likelihood  function  of  |  Y(k  \  ,k2)  \ ,  from  the 
relationship 

Y(k],k2)  =  l(k^k2)-X(k],k2).  (4.3.1. 1) 

Then,  since  the  noise  sequence  C,'(mx,m2)  is  assumed  to  be  white  Gaussian, 
W(k  i  ,k2)  and  Y(k\,k2)  follow  the  Rayleigh  distribution  as  in  the  following  theorem. 
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Theorem  4.3.1:  The  modulus  of  the  DFT  of  the  noise  sequence,  {\W(k\,k2)\, 
ky=0,l,...fl-l,  k2=0,l,...,  N/2  }  and  {\Y{kx,k2)\,  *1  =0,1 . N-l, 


k2=0, j/V  /2  \j  are  the  white  sequences  with  the  following  Rayleigh  densities. 


Z,  rr  ,  rr  j 

f\w{kl.k2)\(W)  =  <-^eKPpN2  ,W>0 

0  ,  otherwise 


2s\lMY{kuk2)  t  ~s2kl,k2Y2(k\M) 
fmtM,  O')  =  |  ^  exP  t  ^2  J  ,  Y(k ,  ,*2)  >  o 

0  ,  otherwise 


where 


d_ 

2nn\  T  2ji/J2  2 

S*„A2  =  I  2 COS(-  )  -  26-05  (0!  I  1 2 cos(-~N  )  -  2co.9(02  I 


{Proof)  Let  W'(ki  ,k2)  be  the  DFT  of  a  white  noise  sequence  C  (mi  >m2)  with  a  normal 
distribution,  that  is. 


£'(m,,m2)  ~ /V(0,p)  ,for  irq ,m2  =  1,2,...,N 


(4.3. 1.2) 


Then,  W'(k!  ,k2)  will  follow  another  normal  distribution  as  follows  [26]. 


W'(k,  ,k2)  -  /V(N-E[C'(m,  ,m2)1,pN2) 


(4.3.1.3.a) 


/V(0,pN2) 


(4.3. 1. 3. b) 


and 
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pN2 

Re(W'(k1,k2)),  Im(W'(ki  ,k2))  -  N (0,  — )  (4.3. 1.4) 

where  Re(W')  and  Im(W')  are  the  real  and  the  imaginary  pans  of  W'(k],k2), 
respectively. 

Define  new  coordinate  systems,  (1],12)  and  (n1(n2),  which  can  be  obtained  from 
the  coordinate  transformations  (4.2. 3.2.5),  (4.2.3.2.10).  Thus,  based  on  the  new 
coordinate  system,  a  new  noise  sequence  £(li,l2)  vvhich  has  the  corresponding  DFF, 
W(kj,k2),  can  be  defined.  As  defined  before,  since  the  density  function  of 
Re(W'(ki ,k2))  or  Im(W'(k! ,k2))  follows  a  normal  distribution  and  the  coordinate 
transformations  are  the  linear  transformations,  the  density  function  of 
^(li ,l2),Re(W(k! ,k2))  will  follow  the  same  normal  distributions  to  the  ones  of 
^'(m1,m2),  Re(W'(k[  ,k2)).  That  is, 

C(li,l2)~(V(0,p)  ,for  1,  ,12  =  1,2,...,N  (4.3.1.5) 

and 

Re(W(ni,n2)),  hn(VJ(nun2))~N(0,  £^1)  (4.3.1.6) 

Now,  consider  the  density  function  of  |W(n!,n2)|.  Since  W(n!,n2)  is  a  complex, 
defining  W  be  W(nt  ,n2), 

|W(n,,n2)|  =  VRe2(W)  +  Im2(W)  (4.3.1.7) 

Define  <1)  =  tan-1  (—  ^— ^-).  Then,  we  can  have  the  following  relations. 

Y  Re(W)  b 

Rc(W)  =  j  W(n,  ,n2)  I  cos<|>  (4.3.1.8.a) 

lm(W)=  |W(n1,n2)|  sin<j>  ,for|W|  >0  (4.3. 1. 8. b) 

Thus,  the  joint  probability  density  of  |W'|  and  $  can  be  obtained  from  the  Jacobian 
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and  the  joint  probability  density  of  Re(W')  and  Im(W')  as  follows. 

/|W|  <*>(l W|,<J))  =  -j— ^T^-/Re(W)Im(w)(Re(W),Im(W))  ,for|W|  >0 

=  I W  |/rC(W)(  I W  |  cos<j>)  /!m(w)(  1 W  |  sin<t>) 
-(Re2(WHIm2(W)) 

_  iwi  c - - 

7tpN2 


(4.3.1.9.a) 

(4.3.1.9.b) 

(4.3.1.9.C) 


where  Jacobian  = 


0|W| 

0Re(W) 

0<j) 

0Re(W) 


0 1 W  | 

0Im(W) 

0<j) 

0Im(W') 


(4.3.1.10.a) 


1 

|W| 


(4.3-1.10.b) 


Therefore,  since  |  W|  and  $  are  independent,  and  the  density  function  of  (|),  /,*,((])),  is 
1 

2n’ 


/|W|o(|W|,phi)=/|W|(|W|)-/0(4>) 


(4.3.1.11) 


and 


/|W|(W)  = 


2W  .  -W2  , 

- n-exp( - r) 

pN2  pN2  ,  W>0 
0  ,  otherwise 


(4.3.1.12) 


Now,  set 


c  A 

27tni  T  27tn2  2 

sk,.k2  =  1 2cos(  —  )  -  2coscoi  |  |2cos(-  —  )  -  2cosco2 1  (4.3.1.13) 

Then,  from  the  definition  of  nj,  (2.3.2.10),  we  can  have 
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|  Y(kj  ,k2)  |  =skl>kl-1  !  W(k!,k2)|  (4.3.1.14) 

=  skl,k2‘1  I  W'(n,  ,n2)  |  (4.3.1.15) 

Thus,  the  density  function  of  |Y(k1,k2)|  can  be  represented  by  following 
equations. 

f|Y(kj,k2)|  00  =  skltk2  f|W(k,.k2)|(sk,.k2Y(kl.k2))  (4.3.1.16) 


2s2k„k2Y(k1,k2)  -s2kl.k2Y2(k1,k2) 

- x - exp  {  - 2 - ) 

pN2  pN2 


,  Y(k,  ,k2)  >  0 
,  otherwise 


From  these  probabilistic  properties,  the  estimation  of  parameters  can  be  done  by 
a  hybrid  method  of  Least  Square  and  Maximum  Likelihood  estimations,  which  was 
suggested  by  Eom  [67].  For  LS  estimation,  if  the  values  of  the  illumination  direction 
cjl,  Tl  and  the  normalization  factor  e  are  given,  and  parameters  tOj,  co2,  o,  x  are  set, 
then  0  =  (c,d,a)T  can  be  estimated  by  minimizing  the  following  cost  function. 

N-l  H 

J(e,<B1,co2)=  21  E  (log  |  I(kj  ,k2)-N  e[sinocosxcosTLsinaL 

ki=0  k2=0 
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d  2jt[(cosa-l)sinxcosxkj+(cos2x+cosasin2x)k2]  7 

+  —  l°g  1 2cos - — - 2coso)2  I  +  a) 

2  Ncosa 


(4.3.1.17) 


=  £  (log  I  I(k!  ,k2)  -  N2e[sinacosxcosxLsinaL+sinasinxsinxLsinaL 

k,=o  k2=0 


+cosCTcosaL]8(ki,k2)|  - 0TQ(ki,k2))2  (4.3.1.18) 

where  a  =  -E|log  |  W(k,  ,k2)  |  ],  0  =  (c,d,a)T. 


Here, 


Q(ki,k2)  = 


-y-  log  1 2cos 
—  log(2cos 


27t[(sin2x+cosacos2x)k1-f(cosa-l)sinxcosxk2] 

Ncosa 

2jt[(cosa~l)sinxcosxki+(cos2x+cosasin2x)k2] 

Ncosa 

-1 


-  2cost0!  I 

-  2cos&)2  I 


(4.3.1.19) 


Thus,  the  estimated  values  will  be 


N/2 


N/2 


fc,d,a]T  =  (N£  LSJQ(k1,k2)QT(k1,k2)r1(N2  L2JQ(k1,k2)log|Z(k1,k2) 


k,=0  k2=0 


k,=0  k2=0 


-  N2e[sinacosxcostLsinaL+sinasinxsinTLsinaL+cosacosaL]5(k1 , k2)  I )  (4.3. 1 .20) 

Also,  ML  estimators  of  (olt  002,  a,  x  can  be  calculated  by  maximizing  the  log- 
likelihood  function  L(Y;0,O)i  ,02)  with  a  estimated  from  above. 
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N-l  M 

L(Y;0,co1,co2)  =  £  £  log|I(k,,k2)-N2e[sinacosxcosxLsinaL 


k,=0  k2=0 


N-l  k 


+  sinasintsinXLsinOL  +  cosacosajj5(k,,k2)|  -  £  Z  ^°g(~~~  ~) 

k,=0  k2=0  2 


c  N-i  27t[(sin2x+cosacos2x)k,+(cosa-l)sinxcosxk2] 

+  —  N  £  log  1 2(cos - - - cosco, )  | 

2  Ncoso 


[N/2.  2k[(coso-1  )sinxcosxki+(cos2x+cososin2x)k2] 

+  dN  £  log|2(cos - — - cosco2)| 


Ncosa 


l  n-i  N/2j  2jc[(sin2x+cosGcos2x)kt+(cosG~l)sinxcosxk2] 

— y  Z  Z  1 2 (cos - - — - - cosco,) | 

pN2  k^o  Ncos° 


27r[(cosa-l)sinTcosxkj  +(cos2x+cosasin2x)k2] 

|2(cos - — - cosco2)  | 

Ncosa 


|  I(k,  ,k2)-N2efsinacosxcosXLsinaL+sinasinxsinxj_sinaL+cosacosaL]6(ki  ,k2)  | ' 


(4.3.1.21) 


where  p  is  a  variance  of  £(1,,12)  and  can  be  estimated  by  the  following  equation 
in  the  mean  square  sense. 


(4.3.1.22) 


where  Y  is  Euler’s  constant  (=  0.5772157)  [67j. 


Therefore,  the  estimation  scheme  can  be  summarized  as  follows: 
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In  this  experiment,  a  whole  image  which  contains  the  shade  and  texture  on  a  3-D 
sphere  surface  was  constructed  using  our  proposed  3-D  texture  model,  and  this  was 
compared  with  the  images  which  were  obtained  by  either  applying  the  reflectance  map 
function  or  by  projecting  the  texture  pattern  only.  Since  our  mode!  is  a  composite  of 
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the  Shape-from-shading  and  the  Shape-from-texture  models,  the  constructed  image 
will  look  more  natural  than  the  ones  from  either  the  Shape-from-shading  or  the 
Shape-from-texture  technique.  For  this  experiment,  the  texture  pattern  was  chosen  as 
the  Fractional  Differencing  Periodic  model  (4.2.3. 2.6)  with  the  parameter  values, 
C0i  =  0.2,  g>2  =  0.2,  c  =  0.8,  and  d  =  0.8.  From  the  given  illumination  direction  oL, 
and  the  surface  orientations  of  sphere,  a,  x,  the  3-D  surface  of  a  sphere  covered  by  the 
chosen  texture  pattern  was  synthesized.  Here,  the  illumination  direction  was  chosen  as 
Ol  =  -0.66,  tL  =-0.66,  and  the  surface  orientations  o(i},i2),  T(ij ,i2)  of  the  (i] , i2)th 
patch  was  given  by 


o(ii,i2)=  1 


tan-1 ( 


a//ai.i2)/ai2 
d//(i,,i2)y3i1  li,,i2 


o 


,if//(l!,l2)>63 

.otherwise 


(4.3.2.1) 


and 


x(ii,i2)  = 


cos 


0 


-l 


^0///ai1)2|il,i2  +  0H/B 12)2  Ii„i2  + 


,if  //(lj  ,12)  >  63  (4-3.2.2) 
.otherwise 


where  the  heig!  ♦  function  H{\\ ,12)  is  given  by 

//(li,l2)=  152  -  lj2  -  122  (4.3.2. 3) 

In  this  experiment,  each  3-D  texture  pattern  was  synthesized  on  32  x  32  pixel  sized 
planar  patch,  and  16x16  pixel  sized  patch  was  taken  from  the  center.  The  complete 
image  of  sphere  (512  x  512)  was  obtained  by  adjoining  these  16  x  16  pixel  sized 
patches.  Figure  4.3-a,b  show  the  3-D  shape  of  a  hemisphere  and  the  corresponding 
sphere  image  obtained  by  the  reflectance  map  function  (4.2. 2. 3)  with  the  parameter 
values  Cl  = -0.66,  Tl  = -0.66,  e  =  100  (4.3.2),  and  Figure  4.4  shows  the 
orthographically  projected  texture  image  based  on  the  given  sphere  surface  with 
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equation  (4.2. 3. 2. 6).  Finally,  with  the  composite  model  of  the  projected  texture  image 
and  the  reflectance  map  function  (4.3.1),  an  image  of  3-D  texture  on  the  surface  of  a 
sphere  was  synthesized  (Figure  4.5). 

Here,  Figure  4.3  and  Figure  4.4  illustrate  that  neither  the  Shape-from-shading  nor 
the  Shape-from- texture  technique  is  suitable  for  representing  the  natural  scene.  Figure 
4.3  does  not  contain  any  detail  information  about  the  surface  pattern,  while  Figure  4.4 
does  not  contain  any  shade.  Thus  it  is  difficult  to  see  the  sphere  in  either  Figure  4.3  or 
Figure  4.4.  Figure  4.5,  however,  illustrates  the  sphere  image  with  texture  pattern  on 
its  surface  clearly,  thus  shows  our  model’s  capability  to  represent  the  3-D  textural 
surface. 

4.3.2.2.  Parameter  Estimation 

From  this  experiment,  we  want  to  show  how  accurate  the  estimated  results  are. 
Since  each  small  patch  is  assumed  to  be  a  tilted  and  slanted  texture  plane  and  the 
whole  image  is  obtained  by  adjoining  these  patches,  our  proposed  3-D  texture  model 
will  fit  each  small  image  patch  and  the  surface  orientation  parameter  will  be  estimated 
based  on  each  patch.  Thus,  in  this  experiment,  we  will  consider  single  patches  of 
texture  patterns  which  represent  the  tilted  and  slanted  texture  planes.  From  these 
patches  the  parameter  values  of  the  model  (4.3.1)  will  be  estimated  by  the  proposed 
estimation  scheme,  and  compared  with  the  true  values.  For  this  experiment,  three 
different  2-D  texture  patterns  sized  64  x  64  were  generated  by  equation  (4.3.1)  with 
the  different  values  of  parameters  G)j ,  0)2,  c,  and  d.  Then,  the  projected  images  of  the 
slanted  and  tilted  texture  planes  were  synthesized  with  the  different  values  of  the 
surface  orientation  a,  x  in  equation  (4.3.4),  and  zero  mean  white  Gaussian  noise  with 
variance  10  was  added  to  each  image.  The  values  of  the  illumination  direction  Cl,  Xl 
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(a) 


Figure  4.3:  Sphere  images:  (a)  Height  function  of  a  sphere  obtained  by 
equation  (4. 3. 2. 3)  (b)  Image  obtained  by  the  reflectance  map  function 
(4. 2.2. 3)  with  £=100  (4.3.2),  G|.  =  -0.66,  xL  =  -0.66 
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Figure  4.4:  Image  obtained  after  projecting  texture  image  orthographically 
to  the  sphere  surface.  (Background  texture  pattern  is  generated  by 
Fractional  Differencing  Periodic  model  (2. 3. 2. 6)  with  o)i=0.2,  co2=0.2, 
c=0.8,  d=0.8) 


'JO 


Figure  4.5:  Image  of  3-D  texture  on  the  surface  of  a  sphere  (Image  is 
generated  by  the  composite  model  (4.3.1)  with  Gl  =  -0.66,  Tl  =  -0.66) 
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and  the  normalization  factor  e  was  given  by  — ,  — ,  and  100,  respectively.  Thus,  since 

8  4 

the  intensity  value  of  a  plane  surface  will  be  constant  all  over  the  patch,  the 
deterministic  part  of  equation  (4.3.4)  was  considered  as  a  constant,  and  added  to  the 
random  field  part.  For  the  estimation,  the  hybrid  method  of  Least  Square  and 
Maximum  Likelinood  estimation,  which  was  discussed  in  previous  section,  was  used 
to  estimate  parameters. 

Figure  4.6-a  shows  one  of  the  2-D  texture  pattern  which  was  generated  by 
equations  (4.2.3.2.6)-(4.2.3.2.12)  with  the  values  of  CO],  CO2,  c,  d  as  0.2,  0.2,  0.8,  0.8, 
respectively,  and  Figure  4.6-b  shows  the  projected  image  of  the  version  of  Figure  4.6- 

a  tilted  and  slanted  by  ~  and  respectively.  As  can  be  seen  from  Figure  4.6-a  and 


Figure  4.6-b,  due  to  the  slanting,  the  distance  between  the  dark  spots  along  the  normal 
direction  of  the  tilted  axis  (the  tilted  axis  is  the  North-west  to  South-east  diagonal) 
reduces  considerably  (thus  producing  almost  a  continuous  dark  band  along  that 
direction  in  Figure  4.6-b).  Here,  the  2-D  texture  model  does  not  fit  the  pattern  in 
Figure  4.6-b,  since  this  pattern  has  a  non-isotropic  random  texture  distribution  due  to 
the  tilt  and  slant.  Note  that  this  synthesis  technique  does  not  require  an  interpolation 
after  projecting.  Because  the  white  Gaussian  random  noise  will  still  be  white  Gaussian 
noise  after  projecting,  the  random  noise  as  input  data  can  be  generated  for  the 
projected  model  directly.  From  these  three  different  3-D  texture  patches,  we  have 
estimates  close  to  the  tme  values.  [Table  4.1]  Also  shown  in  the  parentheses  are 
absolute  deviations  from  the  true  values.  Therefore,  we  can  say  that  the  height 
function  and  the  texture  pattern  of  whole  image  can  be  estimated  properly  from  this 
local  patch  estimation  process. 
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(a) 


(b) 


Figure  4.6:  2-D  texture  images:  (a)  Image  obtained  from  (4.2. 3.2.6)  with 
w^O.2,  <*>2=0.2, c=0.8,d=0.8  (b)  Projected  image  of  the  tilted  and  slanted 
version  of  Figure  4.6-a  by  -7t/4  and  7t/8,  respectively 
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Table  4.1:  True  parameter  values  and  the  estimated  parameter  values  of  the 
projected  Fractional  Differencing  Periodic  model  (4.2.3.2.6):  The  random  noise 
sequence  as  input  data  for  each  synthesized  texture  patch  of  size  64x64  was 
generated  from  white  Gaussian  noise  with  zero  mean  and  variance  10.  Estimated 
values  were  obtained  from  these  synthesized  images  by  equations  (4.3. 1.8)- 
(4.3.1.13). 


Patches 

True  Values 

CO, 

o? 

c 

d 

G 

X 

Patch  1 

0.2 

0.2 

0.8 

0.8 

“0.4 

0:78 

Patch  2 

0.2 

0.2 

0.4 

0.8 

0.78 

0.78 

Patch  3 

0.2 

0.4 

0.6 

0.8 

0.2 

0.6 

4.3.2.3.  Tree  Image 
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A  512x512  real  image  of  a  part  of  a  tree,  which  contains  shade  and  texture  on 
surface  was  considered  (Figure  4.7-a).  From  the  whole  image,  the  estimation  values  of 
the  illumination  direction  Gl,  Tl  were  determined  to  be  Gl  =  -1.3384,  =-0  0426, 

by  applying  Lee’s  method.  Then,  for  this  experiment,  several  local  patches  of  size 
32x32  were  taken  from  the  tree  surface,  and  the  parameter  values  were  estimated  by 
applying  our  projected  Fractional  Differencing  Periodic  model  (4.3.1)  and  estimation 
scheme  (4.3. 1.7)-(4.3. 1.13).  Average  estimated  values  of  the  parameters  were 
obtained  as  coj  =  0.486,  0)2  =  0.053,  c  =  1.394,  d  =  0.762.  A  Synthesized  tree  surface 
image  with  the  given  illumination  direction  was  constructed  by  projecting  this  texture 
pattern  on  a  cylinder  surface  (Figure  4.7-b).  This  image  looks  very  similar  to  the 
original  tree  image,  and  shows  the  ability  of  our  model  to  represent  a  non-stationary 
texture  pattern  such  as  a  tree  bark  on  tree  surface. 

4.4.  Conclusions 

In  this  chapter,  a  composite  model  of  Shape-from-shading  and  Shape-from- 
texture  based  on  a  2-D  orthographically  projected  Fractional  Differencing  Periodic 
model  was  developed  to  represent  a  3-D  surface  image  which  contains  information 
about  both  radiance  and  texture.  This  composite  model  has  several  advantages  over 
the  conventional  approaches.  First,  as  compared  to  the  Shape-from-shading 
techniques,  this  model  always  gives  unique  and  more  accurate  solutions  for  the 
surface  orientation  parameter  values,  because  of  the  additional  constraint  from  the 
texture  function  part.  Also,  by  using  this  analysis,  the  integrability  problem  which 
might  occur  in  spatial  domain  analysis  can  be  avoided,  because  only  one  inverse 
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(a) 


(b) 


Figure  4.7:  Tree  images:  (a)  A  512x512  original  image  of  a  part  of  a  tree,  (b)  A 
synthesized  tree  image:  Local  patch  size  is  32x32,  and  each  3-D  texture  pattern 
was  synthesized  by  the  composite  model  (4.3.1)  with  the  illumination  direction  (?l 
=  -1.3384,  tl  =  -0.0426  and  coj  =  0.486,  co2  =  0.053,  c  =  1.394,  d  =  0.762  for  the 
random  part.  16x16  pixel  sized  patch  was  taken  from  the  center  of  it.  The 
complete  image  of  cylinder  (512x512)  was  obtained  by  adjoining  these  16x16 
pixel  sized  patches. 


Fourier  transform  needs  to  be  taken  at  the  end  of  procedure  to  get  the  whole  image. 
Second,  as  compared  to  the  Shape-from-texture  techniques,  the  Fractional 
Differencing  model  has  the  property  of  being  flexible  enough  to  explain  both  the 
long-term  and  the  short-term  correlation  structure  of  the  texture  pattern,  thus,  it  has  a 
superior  ability  to  model  different  textures  encountered  in  practice.  The 
orthographical  projection  adds  the  additional  flexibility  to  represent  the  3-D  rotated 
texture  due  to  the  slant  and  the  tilt  of  a  surface  normal  plane.  The  estimation  scheme 
for  the  parameters  was  based  on  the  hybrid  method  of  least  square  and  maximum 
likelihood  estimations. 
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CHAPTER  5 
CLASSIFICATION  OF 
3-D  ROTATED  TEXTURES 


5.1.  Introduction 

Texture  classification  has  been  the  focus  of  interest  to  many  researchers 
[53,92,126,216J.  The  classification  problem  can  be  stated  as  allocation  of  an 
observed  texture  image  data  to  the  one  of  the  pre-defined  texture  classes.  These 
texture  classes  can  be  described  by  texture  features,  and  then  texture  features  can  be 
the  parameters  in  stochastic  [53,92,126]  or  structural  models  [216].  Thus,  the  key 
step  in  classification  process  is  the  choice  of  a  set  of  features  which  can  reduce  the 
dimension  of  the  image  data  to  a  computationally  reasonable  amount  of  data. 
Preferable  featuies  should  be  those  that  are  simple  and  easy  to  extract  from  the  given 
data  while  preserving  the  classifying  information  present  in  it. 

Most  of  classification  schemes  which  have  been  suggested  up  to  date  are  under 
the  assumption  that  the  test  sample  data  possesses  the  same  surface  orientation  as  the 
training  sample  data.  Thus,  if  the  orientation  of  test  image  is  different  from  the 
training  sample  data,  for  example,  in  case  of  a  rotated  image,  the  classification 
performs  poorly.  This  reduces  the  flexibility  of  those  classification  schemes. 
However,  most  of  natural  texture  images  which  we  encounter  in  practice  represent  the 
texture  of  3-D  surfaces.  Thus,  the  observed  image  is  a  projected  surface  image  onto 
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the  2-D  image  plane  with  a  2-D  rotated,  or  3-D  slanted  and  tilted  texture  pattern. 
Therefore,  it  is  desirable  that  the  classification  scheme  have  the  flexibility  to  classify 
even  rotated  or  scaled  texture  to  the  original  class  of  it.  This  is  a  good  indication  why 
it  is  so  important  to  have  the  rotational  and  scaling  invariant  features  in  our  model. 
Up  to  date,  several  approaches  have  been  reported  for  the  rotational  invariant 
classification  scheme  [53, 126].  However,  those  approaches  have  their  own  limitations 
of  classification.  Khotanzad’s  method  [126]  is  based  on  the  four  different  discrete 
directions  instead  of  an  arbitrary  angle  of  the  rotational  direction,  thus,  with  this 
classification  scheme,  the  rotated  texture  pattern  with  an  arbitrary  angle  can  not  be 
classified.  Cohen’s  model  [53]  uses  a  coordinate  transformation  for  the  rotation  and 
the  slant  of  texture  surface,  and  a  ML  estimation  method  is  adapted  to  extract  the 
texture  pattern  parameters  in  an  AR  model  and  the  surface  orientation  parameters 
simultaneously.  However,  because  of  the  complexity  of  the  estimation  scheme  and  the 
limitation  of  AR  model  to  represent  the  surface,  this  approach  may  not  be  practical.  In 
this  chapter,  a  multi-level  classification  method  which  can  handle  arbitrary  3-D 
rotated  samples  of  textures  is  developed  based  on  fractional  differencing  models  with 
a  fractal  scaling  parameter.  Since  the  fractal  scale  is  known  to  be  a  rotational  and 
scaling  invariant  parameter,  the  accuracy  of  classification  will  not  be  affected  by  3-D 
rotation  of  the  test  texture,  by  using  these  models.  In  first  level  of  classification,  the 
textures  are  classified  by  the  first-order  Fractional  Differencing  model  with  a  fractal 
scale  parameter,  and  in  second  level,  classification  is  completed  with  the  additional 
frequency  parameters  of  the  second-order  Fractional  Differencing  periodic  model. 
This  multi-level  classification  scheme  has  at  least  following  advantages  over  the 
conventional  approaches  [46,53, 126].  First,  since  the  fractal  scale  parameter  of  the 
first-order  Fractional  Differencing  model  can  be  estimated  by  a  simple  Least-square 
estimation  method,  the  processing  time  can  be  dramatically  reduced.  Second,  the 
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ambiguity,  which  may  be  caused  by  using  only  one  classification  parameter, 
particularly  when  the  values  of  parameter  are  close  enough  for  the  different  textures, 
can  be  removed  by  considering  the  additional  parameters  in  second-order  model.  Even 
though  the  estimation  scheme  for  these  additional  parameters  is  based  on  ML 
estimation,  the  classification  process  is  much  simpler  compared  with  other  approaches 
based  on  ML  estimation,  e.g.  Cohen’s  [53],  because  only  a  small  number  of  texture 
patterns  in  same  sub-class,  which  is  already  classified  based  on  the  fractal  scale  value 
in  first  level,  need  to  be  classified. 

The  organization  of  this  chapter  is  as  follows.  In  section  5.2,  the  3-D  rotated 
texture  is  defined  by  the  linear  coordinate  transformations.  Section  5.3  introduces  the 
Fractional  Differencing  model  with  a  fractal  scaling  parameter  to  represent  the  rotated 
and  slanted  texture  surface,  and  section  5.4  discussed  its  parameter  estimation  scheme 
based  on  LS  and  ML  estimation  methods.  In  section  5.5,  a  multi-level  3-D  rotational 
invariant  classification  scheme  is  introduced.  Section  5.6  then  discusses  some 
simulation  results  carried  out  to  demonstrate  the  performance  of  the  proposed 
algorithms,  followed  by  section  5.7  which  concludes  this  chapter. 

5.2,  3-D  Rotated  Texture 

To  represent  a  3-D  rotated  texture,  we  need  two  different  sets  of  coordinate 
transformations.  First  one  is  the  2-D  rotational  coordinate  transformation  and  second 
one  is  the  orthographical  projection  coordinate  transformation.  The  2-D  rotational 
coordinate  transformation  is  needed  to  represent  a  texture  image  rotated  on  the  2-D 
image  plane,  and  the  orthographical  coordinate  transformation  is  needed  to  represent 
an  orthographically  projected  texture  surface  onto  the  image  plane  due  to  the  slant  and 
tilt  of  texture  surface.  Here,  notice  that  for  the  image  whose  texture  pattern  is 
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distributed  isotropically  over  the  surface  without  any  directional  trend,  the  2-D 
rotational  coordinate  transformation  is  not  required,  because  the  2-D  rotation  will  not 
affect  the  pattern  of  the  texture. 

The  coordinate  system  (/' \,l'i)  which  is  rotated  on  the  image  plane  through  an 
angle  0  clockwise  can  be  obtained  from  the  original  coordinate  system  (/ 1  ,/2)  by  the 
following  simple  coordinate  transformation. 

(5.2.1) 


[''ll 

COS0 

-sin9 

r  *1 

/ 1 

I'l 

sinQ 

COS0 

h 

And  a  texture  surface  plane  slanted  and  tilted  by  a  and  t  can  be  represented  by  the 
orthographical  projection.  First,  take  the  mj-mj  coordinate  system  over  the  surface 
plane.  Put  a  line  passing  through  the  origin,  and  let  T  be  the  angle  made  from  the  m  j  - 
axis.  Rotate  the  plane  around  the  line  by  angle  a  and  project  the  rotated  plane 
orthograph  ically  onto  the  image  plane.  Thus,  a  new  coordinate  system  from  the 
viewing  direction,  l\-l2,  can  be  obtained  from  the  following  two  coordinate 
transformations. 


V 

cost  -sinT 

m\ 

111 

sint  cost 

m2 

and 


V 

1  0 

V 

l7 

0  cos 0 

- 

111 

(5.2.2) 


(5.2.3) 


Hence,  as  is  well  known  [110),  the  coordinate  transformation  of  the  orthographic 
projection  between  the  m  i  -m2  system  and  / 1  -l 2  system  can  be  given  as  follows. 


“ 

* 

- 

■ 

-1 

- 

1 1 

cost  -sinT 

1  0 

cost  -sinT 

m, 

1 2 

sinT  cost 

0  cosa 

sinT  cost 

m2 

(5.2.4a) 


104 


cos2x+cosasin2x  ( l-coso)sinxcosx 
(l-coso)sinxcosx  sin2x+cosacos2x 


m\ 

m2 


(5.2.4b) 


Therefore,  a  texture  plane  coordinate  system  (l'\,l"2)  obtained  after  rotating, 
slanting  and  tilting  with  arbitrary  angles  can  be  represented  by  combining  both  2-D 
rotational  coordinate  transformation  (5.2.1)  and  the  orthographical  projection 
coordinate  transformation  (5.2.4)  as  follows. 


'/"t ' 

cos0 

-sin0 

cos2x-i-coscsin2x 

(l-cosa)sinxcosx 

m  1 

/"  2 

sin0 

COS0 

(l-cosa)sinxcosx 

sin2x+cosocos2x 

m2 

(5.2.5a) 


/  II 
~  III  IV 

where 

I  =  cos0(cos2x  +  cosasin2x)  -  sin0(l-cosa)sinxcosx 


(5.2.5b) 


(5.2.6a) 


II  =  cos0(  l-coso)sinxcosx  -  sin0(sin2x  +  cosacos2x)  (5.2.6b) 

III  =  sin0(cos2x  +  cosasin2x)  +  cos0(l-cosa)sinxcosx  (5.2.6c) 

IV  =  sin0(  l-coso)sinxcosx  +  cos0(sin2x  +  cosacos2x)  (5.2.6d) 


One  grid  pattern  image  (Figure  5.1 -a)  was  considered  to  demonstrate  these  coordinate 
transformations.  The  coordinate  transformation  was  taken  to  this  image  with  0  =  7t/8  , 
a  =  7t/4,  and  x  =  -rc/8.  Figure  5.1-b,  Figure  5.1-c,  and  Figure  5.1-d  depict  a  2-D 
rotational  coordinate  transformation,  a  orthographical  projection  coordinate 
transformation,  and  both  coordinate  transformations,  respectively. 
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Figure  5.1:  2-D  grid  pattern  images:  (a)  Original  image,  (b)  Image  obtained  after 
rotating  the  image(Figure  5.1 -a)  on  the  image  plane  with  0  =  7t/8  in  (5.2.1),  (c) 
Image  obtained  after  projecting  the  image(Figure  5.1-a)  orthographically,  with 
0=71/4  and  x  =  -7t/8  in  (5.2.4),  (d)  Image  obtained  after  projecting  the  already 
rotated  image(Figure  5.1-b)  orthographically,  with  o  =  7t/4  and  T  =  — rc/8  in  (5.2.5). 
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5.3.  Fractional  Differencing  Model  With  One  Fractal  Scaling  Parameter  c. 

All  features  considered  in  this  chapter  for  the  classification  purpose  can  be 
obtained  by  fitting  a  two-dimensional  parametric  random  field  models  to  the  given 
texture.  There  are  various  kinds  of  random  field  models  suggested  up  to  date 
[4,41,56,85,116],  and  each  random  field  model  has  its  own  advantage  and 
disadvantage  based  on  the  purpose  of  the  process.  Among  the  various  random  field 
models,  the  fractional  differencing  models  with  one  fractal  scaling  parameter  are 
chosen  in  this  study  because  of  its  following  properties.  First,  the  fractal  scaling 
parameter  c  is  known  to  be  a  rotational  and  scaling  invariant  parameter,  thus,  with  this 
parameter  we  have  the  flexibility  to  handle  the  rotated,  slanted,  and  tilted  texture 
surface.  Second,  unlike  the  other  fractal  model  based  on  the  fractional  Brownian 
process,  this  model  has  a  simple  estimation  scheme,  such  as  least  mean  square 
estimation  for  the  parameter  c.  Third,  with  the  second  order  periodic  model,  this 
model  has  good  performance  in  texture  synthesis  and  its  ability  to  simultaneously 
represent  the  coarseness  and  pattern  of  texture  surface  with  the  fractal  scaling 
parameter  c  and  directional  frequency  parameters  coi,  0)2,  respectively.  Thus, 
comparing  these  parameter  values,  we  can  classify  the  texture  patterns  properly  even 
though  some  texture  patterns  share^  the  same  value  of  one  of  those  parameters. 
Typical  first-order  and  second-order  fractional  differencing  models  with  one  fractal 
scaling  parameter  will  be  as  follows,  respectively. 

_  C 

y(mi,m2)  =  [(1-Z/-1)(1-Z2_1)]  2t,(mum2)  (5.3.1) 


and 


y(mi,m2)  =  [(\-2cos<olzi  l+zx  2)(\-2cos(02z2  !+z2  2)]  2  ^{mum2)  (5.3.2) 

form!, m2  =  0,1,...,N-1. 

The  corresponding  DFTs  of  these  functions  are 

Y(kl,k2)  =  [(\-e  '  "")(l-e  jlU  N  )]  YW(kuk2),  (5.3.3) 

and 

-;2n— 

Y(k\,k2)  =  [(l-2c<?.s(0ie  N +e  N) 

^*2  ^2  C 

-;2n— -  ~J4n— 

•  (l-2co50)2e  *+<?  *  )]  2  WC*!,^),  (5.3.4) 

where  z,  is  the  delay  operator  associated  with  m,,  £(mi,m2)  is  an  i.i.d.  Gaussian 
sequence,  and  VK  (A:  t  ,fc2)  is  the  corresponding  DFT. 

Note  that  this  random  function  y(mj,m2)  is  defined  on  the  surface  normal  plane 
because  this  model  is  based  on  a  2-D  texture  model.  Thus,  it  will  fit  the  texture 
pattern  which  is  isotropically  distributed  over  the  plane  and  this  isotropically 
distributed  texture  pattern  can  be  obtained  by  viewing  from  the  surface  normal 
direction.  Then,  the  model  for  the  3-D  rotated  texture  pattern  can  be  obtained  by 
applying  the  2-D  rotational  and  the  orthographical  projection  coordinate 
transformations. 

5.3.1.  Rotated  And  Projected  Fractional  Differencing  Model 

The  random  function  y,(/i,/2),  defined  on  the  viewing  direction  image  plane, 
can  be  obtained  by  projecting  y(mx,m2)  onto  the  image  plane  orthographically.  On 
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the  other  hand,  the  rotated  random  function  y"(l\,l2)  can  be  obtained  by  rotating  the 
original  coordinate  axes  through  the  angle  0.  Thus,  if  /i  and  1 2  are  the  axes 
transformed  from  mx  and  m2  based  on  x  and  o,  from  the  coordinate  transformation 
for  the  orthographical  projection  (5.2.4),  we  can  have  the  following  first  or  second- 
order  model. 

c 

/(/1,/2)=[(l-z/r1)(l--?2"1)1  2  C(/t./2)  (5.3.1. 1) 

or 

_  c 

y\l\A2)  =  [(<\-'lcos(iiiz'rx+z\~2)(\-2cos(i}2z'2~l+z'2~2)]  2  C(^1^2)  (5. 3.1.2) 

for/^/2  =0,1,...,N-1, 

where  z\  is  the  delay  operator  associated  with  /,,  and  ^i.l\,l2)  is  an  i.i.d.  Gaussian 
sequence. 

By  the  definition  of  DFT,  W'{k\,k2)  corresponding  to  white  noise  sequence 

C'(/i,/2)is 

N- 1  N-l  2ir 

W'(k1,k2)  =  X  X  ^Vi,h)exp[-j~(mlkx+m2k2)]  (5.3.1.3) 

»«[={)  m2=0 

and,  therefore  from  the  coordinate  transformation  (5.2.4), 

W'(k  i  ,k2)  =  X  X  C'(/ 1  ,h)e*P  H  „ ^ ~  •  [((sin2 t+cosgcos2 x)k  x 

m,=0  m2=0  /VCOSO 

+  (cosa-\)sinxcosxk2)l  x+((coso-l)sinxcosxk  i 

+  (cos2t-f-cososin2x)^2)^2]J  (5.3. 1.4) 


Thus,  we  define 
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sin2x  +  cosccos2t  (cos a  -  1  )simcost 


z  1=Z\ 


(5.3.1.5a) 


(cos a  -  l)sinxcosx  cos2t  +  cosgsm2? 


Z  2  =  Zl 


(5.3.1.5b) 


nx  =  — \ — [(sin2T  +  cosacos2t)/i: !  +  (cosa-  1)^/«tcost^2J,  (5.3.1.6a) 


n2  - [(cosa  -  \)sinTcosxk\  +  (cos2 x  +  cosasirr x)A:2]  (5.3.1.6b) 

cosa 


Then,  from  the  orthographical  projection,  we  can  set 


WXkuk2)  =  W(nu  n2). 


(5.3.1.7) 


Therefore,  with  these  relations,  we  can  have  the  projected  version  of  Y(kx,k2), 
which  is  the  DFT  of  y\l\,l2)  as  follows. 


Y'{kx,k2)  =  [^-e  N  )d-e 


-jln^j-  ~ 

N  )]  2  W'{ki,k2) 


(5.3.1.8) 


-jin——-  -j4n- — 

Y’{kx,k2)  =  [{\-?.cosuxe  N  +e  N) 


•(\-2cosco2e 


"2  .,  n2  c 

]2k-  -,4n_  -T 


)]  W\ki,k2)  (5.3.1.9) 


On  the  other  hand,  the  random  function  /'(/], /2),  rotated  version  of y(l \,l2),  can  be 
obtained  by  the  coordinate  transformation  (5.2.1)  and  setting y"(l\, l2)  =  y(l'\,l'i)- 


/'(/i,/2)  =  [n-z"r')(i-zV1)i  2  ?"(/,, h) 


(5.3.1.10) 


or 
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_  c 

y'ViJ2)  =  [0-2cos(olz',rl+z'\-2)(\-2cos(^z'Yl+z'Y2)]TC,,Vul2X5.3A.n) 

where  z",  is  the  delay  operator  in  direction  of  /',  and  £"(/ 1 ,l2)  is  an  i.i.d  Gaussian 
noise  sequence. 

Like  earlier,  from  the  definition  of  DFT,  we  now  define 


_  _  cos9_  sinG 
Z  \  -  Z\  z2 


(5.3.1.12a) 


and 


_  _  -sinG  _  cosG 
Z  2  -  Z\  2  2 


n  \  =  cos0&i+sin0&2 


(5.3.1.12b) 


(5.3.1.13a) 


n'2  =-sin0/:1+cos0/:2  (5.3.1.13b) 

Thus,  the  rotated  version  of  Y (k  j  ,£2).  which  is  DFT  of  y"(l  1  ,/2)  is 

n  1  ti  2  c 

->271-1  ->271— ~  ~ 

Y'Xkuk2)  =  [(\-e  N  )(l-e  N  )]  2  W'Xk.M)  (5.3.1.14) 

or 

n\  n\ 

-j2n— 

Y"(ki,  k2)  =  [{\-2cosmxe  N +e  N) 

„  n2  ..  n’l  C 

-j  271——  —  — — 

•{\-2cosvy2e  N+e  N  )]  2  W"(k}  ,k2)  (5.3.1.15) 

where  W"(k  1  ,k2)  is  the  corresponding  DFT  of  £"(/ 1  ,l2) 

Finally,  the  above  two  coordinate  transformations,  the  2-D  rotational  and  the 
orthographical  projection  coordinate  transformations,  can  be  combined  to  get  the 
model  for  the  texture  pattern  obtained  after  rotating,  and  slanting  the  plane  with  a  tilt 
angle.  Thus,  the  random  function  y'"(l\J2)’  which  represents  the  texture  pattern 
obtained  after  rotating  and  slanting,  can  be  represented  by 
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c 

/"(/i,/2)=[a-z,v1)(i-z/v1)]  Tr(/1,/2) 


(5.3.1.16) 


or 


y"\l i,h)~  [(l-2coA'(D1z"'/_1+z,/,i_2)(l_2co5a)2z",2_1+z,"2-2)l  2^"(lxJ2) 

(5.3.1.17) 

where  z'"  is  the  delay  operator  in  direction  of  /"  and  £"'(/ 1,(2)  is  an  i.i.d.  Gaussian 
noise  sequence. 


And  finally  we  define 


7"'  _  7  Gi  7  G2 

2  1  ~Z1  2  2 

(5.3.1.18a) 

///  "/  21  22 

Z  2=*1  lz2 

(5.3.1.18b) 

where 

T\\  = 

— - — [cos0(sw2x  +  cosocos2!)  -  sin0(co5O  -  1)smtcost], 

COSO 

(5.3.1.19a) 

T 12  = 

1  ■)  9 

- [sin0(sm  T  +  cosocos  x)  +  cos0(coso-  \)sinxcosx], 

COSO 

(5.3.1.19b) 

r21  = 

1  ,  0  'y 

- [cos0(coso  -  ljsmtcost  -  sin0(co^  x  +  cososin  T)], 

(5.3.1.19c) 

COSO 

T  22  ~ 

1  '7  0 

- [sin0(cos 0  -  l)smxcosT  +  cos0(cos  x  +  cososin  x)l 

(5.3.1. 19d) 

COSO 

and 

n"  1  =  Tn/cj  +T  ]2k2. 

(5.3.1.20a) 

n"2  =  T2ikx  +  T22k2 

(5.3.1.20b) 

Thus,  we  get  the  rotated  and  projected  version  of  Y(kx,k 2),  which  is  DFT  of 
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As  shown  in  the  previous  section,  the  random  texture  y(mx,m 2)  defined  on  the 
surface  normal  plane  can  be  viewed  as  the  different  random  texture  depending  on  the 
viewing  direction.  And  this  different  looking  random  texture  can  be  represented  by  the 
coordinate  transformations  with  the  surface  orientation  angles  a,  x  and  the  rotation 
angle  0.  Thus,  the  original  function  y(m  1,^2)  is  transformed  to  another  function, 
such  as  y',  y",  y'",  after  rotating  or  slanting  the  texture  plane.  However,  even  after 
rotated  and  slanted,  the  transformed  random  texture  function  shares  the  same  value  of 
the  fractal  scaling  parameter  c  with  the  original  texture  function  y.  This  rotational  and 
scaling  invariant  property  of  fractal  scaling  parameter  c  plays  an  important  role  for 
our  classification  purpose.  In  other  words,  by  estimating  and  comparing  this  rotational 
and  scaling  invariant  fractal  scale  parameter,  we  classify  the  same  but  different 
looking  textures  due  to  the  3-D  rotation  to  the  same  class.  Following  example  (Figure 
5.2)  is  given  to  demonstrate  this  3-D  rotational  invariant  property  of  the  fractal  scaling 
factor  c.  In  this  example,  all  texture  patterns  are  synthesized  with  different  values  of 
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a,  x,  and  0  but  same  values  of  coj,  o>2,  and  c  based  on  the  transformed  second-order 
fractional  differencing  periodic  models  (5.3.2,  5.3.1.2,11,17).  Figure  5.2-a,  b,  c,  d 
depict  the  original  texture  pattern  with  o)i  =  0.2,  0)2  =  0.2,  and  c  =  0.8,  the  rotated 
texture  with  0  =  71/4,  the  orthographically  projected  texture  with  x  =  tt/8  and  0  =  — 7C/4, 
and  the  rotated  and  projected  texture  with  0  =  71/4,  a  =  -7t/4,  and  x  =  7t/8,  respectively. 
From  this  example,  we  can  see  that  all  texture  patterns  are  same  to  human  eyes,  even 
though  those  are  projected  from  the  different  viewing  directions.  However,  this 
similarity  can  not  be  caught  by  the  conventional  2-D  stochastic  model-based 
approaches  such  as  AR  model  [41,125],  facet  model  [179],  etc.,  but  can  be  caught  by 
the  fractal  scaling  parameter  c  in  fractional  differencing  model,  because  all  texture 
patterns  in  Figure  5.2  shares  the  same  value  of  c.  A  detailed  discussion  on  the 
estimation  method  of  parameter  c  and  the  classification  method  with  this  parameter 
will  be  given  in  the  next  two  sections. 

5.4.  Estimation  Of  Parameters 

There  have  been  several  approaches  for  estimating  parameters  in  various  kinds 
of  fractional  differencing  time  series  since  Hosking  introduced  Fractional 
Differencing  model  [101].  For  example,  Granger  and  Joyeux  [89]  approximated  this 
model  by  a  high-order  auto-regressive  process  and  estimated  the  differencing 
parameter  by  comparing  variances  for  each  different  choice  of  it.  Lapsa  [121] 
suggested  a  maximum-likelihood  estimator  in  the  frequency  domain  and  showed  the 
consistency  of  the  estimator.  This  frequency  domain  analysis  was  further  studied  by 
Eom,  and  a  hybrid  method  of  least-squares  and  maximum-likelihood  estimations  was 
recently  proposed  to  estimate  the  fractal  scaling  parameters  and  the  frequency 
parameters,  respectively  [67].  In  this  section,  a  least- squares  or  a  maximum- 
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(c) 


(d) 


Figure  5.2:  (a)  A  original  texture  pattern  image  with  0)]  =0.2,  o>2  =0.2,  and 
c  =  0.8  in  (5.3.2).  (b)  Image  obtained  after  rotating  the  image(Figure  5.2-a)  on  the 
image  plane  with  9  =  7t/4  in  (5. 3. 1.2).  (c)  Image  obtained  after  projecting  the 
image(Figure  5.2-a)  orthographically,  with  a  =  -7t/4  and  x  =  n/S  in  (5.3. 1.7).  (d) 
Image  obtained  after  projecting  the  already  rotated  image(Figure  5.2-b) 
orthographically,  with  o  =  -7t/4  and  t  -  71/8  in  (5.3. 1.11). 


likelihood  estimation  algorithm  based  on  Eom’s  algorithm  will  be  applied  to  estimate 
the  parameters  in  the  first-order  or  the  second-order  fractional  differencing  model. 

5.4.1.  Estimation  Of  Parameter  c  In  The  First-order  Fractional  Differencing 
Model 

The  estimation  of  fractal  scaling  parameter  c  in  the  un-transformed  first-order 
fractional  differencing  model  (5.3.1)  can  be  done  by  a  simple  least-square  estimation 
scheme  in  the  frequency  domain  based  on  a  representation  of  the  logarithm  of  the 
process  which  is  linear  in  the  parameters  as  follows.  By  applying  logarithm  operator 
to  (5.3.1),  we  obtain 

^1  ^2 

log |  F(& i  ,k2) |  =  -y [log 1 1-e  '  1+  log |  \—e~j2n~^  \  ]  +  log |  IV (k\,k2)  | 

=  -y[log  1 2sin(^~)  |  +  log  1 2sin(^~ -)  |  ]  +  log  |  W(k  j  ,k2 )  \ 

(5.4. 1.1) 

=  — yllog  1 2sin(-^-)  |  +  log  1 2sin(-^-)  |  ]  -  a  +  V(k { ,k2) 

,  for kl,k2  =  0,\,...,N-\  (5.4.1. 2) 

where  a  =  -£[log  I  W(ki,k2)\]  and  V(k\,k2)  =  log |  W(k\,k2)\  +  a. 


Then  r)  =  (c,  a)1  can  be  estimated  by  minimizing  the  following  cost  function. 
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J  (T|,CO],0)2) 


AM  N~ 1  r  Tlk i 

I  £  (iogiy^i.^l  +  -r-  [log  1 2sin— —  | 

*,=0  lc2=0  Z  ™ 


Tik2 

+  log  1 2sin - 

6  N 


|  ]  +  a)2 


=  E  £  OoglK^i,^)!  - Tir<2 (^i,^2))2 

*1=0  *2=0 


Here, 


Q(k  i  ,k 2)  - 


-1  nk ,  nk2 

—  (log  1 2sin— — -  |  +  log  1 2sin— — 
IN  N 


-1 


(5.4.1.3) 

(5.4. 1.4) 


(5.4. 1.5) 


Thus,  the  estimated  values  will  be 

ll*\r  =  (X  NlQ^>k2)QT(kuk2))-'(NZ  NXQ(kuk2)log\Y(k1,k2)\). 

k i  =0  *2~0  *,=0*2=0 


(5.4. 1.6) 


5.4.2.  Estimation  Of  Parameters  CO],  002,  0, t,  0  In  The  Projected  Second-Order 
Fractional  Differencing  Periodic  Model 


Consider  the  transformed  fractional  differencing  periodic  model  obtained  after 
rotating  and  projecting  (5.3.1.17)  and  its  DFT  as  follows. 

n",  n ", 

-j2n~— 

Y'"(klyk2)  =  [(l-2coscoic  N  +e  N  ) 
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-jAn- 


+  e 


N 


)]  lW"\kuk2)  (5.4.2. 1) 


where  W'"{kx,k2)  is  the  DFT  of  y"Xl\,h )  in  (5.3.1.17)  and  n "  is  same  as  defined  in 
(5.3.1.20). 

Then,  the  estimation  of  parameters  in  (5.3.1.17)  can  be  done  in  the  frequency  domain, 
modifying  the  techniques  suggested  in  previous  section. 

For  estimation,  all  parameters  can  be  estimated  directly  from  the  given  data 
Y'”{k\,k2)  (5.3.1.22),  if  we  can  obtain  the  likelihood  function  of  | K(A: j ,A: 2) I  •  Then, 
since  the  noise  sequence  ts"'(l\,l2)  is  assumed  to  be  white  Gaussian,  W'"{kx,k2)  and 
Y"Xk\,k2)  follow  the  Rayleigh  distribution  as  in  Theorem  4.3.1.  From  these 
probabilistic  properties,  the  estimation  of  parameters  can  be  done  by  a  hybrid  method 
of  Least  Square  and  Maximum  Likelihood  estimations,  which  was  suggested  by  Eom 
[671.  For  LS  estimation,  if  parameters  (Oj,  to2,  cr,  x,  6  are  set,  then  t\  =  (c,  a)r  can  be 
estimated  by  minimizing  the  following  cost  function  as  described  in  section  5.4.1. 

N- 1  N~\  c  2nn"i 

7  01,(0!  ,(ri2)=  X  X  (l°glr"(*i>*2)l+y  (log|2ow — - 2COS0)!  1 

*,=0  *2=0  Z  /v 


'l.'Kfl" 

+  log  1 2c 'os — n  -  2cojo)2  | )  +  a)2  (5. 4.2.2) 

=  N£  Nj!(log\Y"Xk,,k2)\-r\TQ(kl,k2))2  (5.4.2. 3) 

*,-0  *2=0 

where  a  =  -E\log  \  W'"(k  j  ,k2)  \  ],  T\  =  (c,  a)7. 
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Thus,  the  estimated  values  will  be 
„  „  ^  N- 1  N- 1 


N- 1  /V-l 


[c,ar=(£  z<2(*i,*2)^|}""(fci,*2)i) 

*,=0*2=0  *,=0*j=0 


(5.4.2. 5) 

Also,  ML  estimators  of  (Oj ,  0)2,  CT,  x,  0  can  be  calculated  by  maximizing  the  log- 
likelihood  function  L(Y with  a  estimated  from  above. 

N- 1  N-l  N-l  N- 1  nA/2 

L(K,T1,0)1,C02)=  X  Elog|y'/'(/:1^2)|-  X 


*,=0  *2=0 
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where  p  is  a  variance  of  C"'(/i,/2)  ar>d  can  be  estimated  by  the  following 
equation  in  mean  square  sense. 

7-2 


"1  *  TiL 

p  =  T^exP  ( y-2a- 


N‘ 


6  N‘ 


(5.4.2. 7) 


where  y  is  Euler’s  constant  (=  0.5772157)  [67]. 
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Therefore,  the  estimation  scheme  can  be  summarized  as  follows: 

(Estimation  Algorithm) 

Step  1:  Choose  resonant  frequencies  tOj,  o>2  in  the  range  [0,7t/2],  the  surface 
orientation  parameters  a,  x  and  the  rotational  parameter  0  in  the  range 
[-71/2, 7C/2] . 

Step  2:  With  the  given  values  of  0)!,  £1)2,  a,  x,  and  0,  estimate  c  and  a  by  LS 
estimation  algorithm  (5. 4.2. 5). 

Step  3:  Using  a,  compute  the  estimate  of  the  variance  of  1 , / 2 )»  P»  by  equation 
(5. 4.2.7). 

Step  4:  Using  the  estimates  c  and  p  found  in  Step  2  and  3,  maximize  the  likelihood 
function  given  by  (5.4.2. 6)  with  respect  to  0)1,  (1)2,  a,  t,  and  0. 

A  A  A  A  A 

Step  5:  Using  the  estimates  coi ,  (1)2,  O,  X,  and  0,  repeat  steps  2  to  4  until  the  estimates 
have  no  significant  change  in  successive  iterations. 

5,5.  Multi-level  3-D  Rotational  Invariant  Classification  Scheme 

In  this  section,  a  two-level  hierarchical  classification  structure  is  developed  to 
classify  the  3-D  rotated  texture  patterns.  In  first  level  of  classification  scheme,  a  3-D 
rotational  invariant  feature  c  is  extracted  from  the  input  image  data  based  on  the  first- 
order  fractional  differencing  model  (5.3.1),  and  this  fractal  scale  feature  is  used  to 
classify  the  textures  to  the  classes  whose  members  are  sharing  similar  value  of  fractal 
scaling  parameter  c.  In  second  level,  each  class  is  divided  to  the  final  desired 
subclasses  based  on  two  other  texture  pattern  features,  O)!  and  o>2,  and  in  each  class 
assigned  in  first  level,  the  input  texture  image  is  classified  further  in  more  detail  with 
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the  features  C0[,  0)2  which  is  extracted  based  on  the  rotated  and  projected  fractional 
differencing  periodic  model  (5.3.1.16).  Figure  5.3  shows  this  hierarchical  structured 
classification  scheme.  Notice  that  it  is  possible  for  some  textures  to  belong  to  the 
different  classes  simultaneously,  depending  their  variance  values  of  c.  For  this 
classification  scheme,  the  images  are  separated  into  test  and  training  sets.  The  class  of 
textures  and  the  number  of  classes  in  training  set  is  assumed  to  be  known  a  priori. 

This  multi-level  classification  structure  has  several  advantages  over  the 
conventional  classification  methods.  First,  this  algorithm  can  save  the  processing  time 
to  classify  the  texture  images.  Because  this  algorithm  uses  a  least  mean  square 
estimation  to  get  t’._  \uctal  scaling  feature  c  in  first-level,  the  processing  time  will  be 
reduced  comparing  with  the  method  which  is  relying  on  only  a  maximum  likelihood 
estimation  techniques.  MLE  algorithm  requires  long  computational  time  because  of  its 
iteration  scheme  to  get  the  maximum  value  from  the  likelihood  function.  And  second, 
comparing  with  the  classification  methods  which  have  only  one  feature  parameter  c, 
the  more  accurate  classification  can  be  achieved  by  this  algorithm.  In  other  words, 
even  for  the  case  that  two  texture  images  have  close  values  of  fractal  scale  (roughness 
of  the  surface)  but  those  are  same  patterns  of  texture,  this  algorithm  has  the  ability  to 
classify  the  textures  based  on  the  estimated  values  of  parameters  (Oi,  (1)2  in  second 
level  of  this  classification  structure.  A  more  detailed  discussion  on  each  level  of 
classification  scheme  will  be  given  in  following  sections. 

5.5.1.  The  First-Level  Of  Classification 

In  this  level,  the  different  3-D  rotated  texture  images  are  classified  into  the  M 
different  classes  depending  on  their  estimated  values  of  the  fractal  scale.  As  discussed 
in  section  5.3.2,  this  fractal  scaling  feature  is  a  rotational  and  scaling  invariant  feature. 
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and  represents  the  roughness  of  the  texture  plane  surface.  In  other  words,  if  we 
consider  a  first-order  fractional  differencing  model,  the  transformed  image  function 
y"Vuh)  obtained  after  rotating  and  projecting  (5.3.1.16)  will  have  same  fractal  scale 
as  the  original  function  v  ( m  \  ,ni2 )  (5.3.1).  Therefore,  if  we  want  to  estimate  the  fractal 
scale  c  only,  we  do  not  need  to  apply  the  transformed  function  but  only  need  to  apply 
the  original  function  which  does  not  contain  the  rotation  parameter  8,  or  the  surface 
orientation  parameters  x,  o  in  it.  Therefore,  the  estimation  of  c  parameter  can  be  done 
by  a  simple  least-square  estimate  method  described  in  section  5.4.1. 

Actual  classification  is  achieved  by  applying  a  distance  classifier  d(c,i),  which 
measures  a  weighted  distance  between  the  extracted  feature  of  test  image  denoted  by 
c  and  the  mean  feature  of  each  of  M  classes.  Then  the  texture  is  classified  to  class  A, 
for  which  such  a  distance  is  minimum.  That  is, 

i*  =  minimum  d(c,i),  i-  1 . M  (5.5.1. 1) 

( 

where 

[c  -  q] 

d(c,i)  —  —jy -  (5.5.1.2) 

Z  [<*c2]0) 

j  =  1 

and  c,  and  [ac2](l)  correspond  to  the  sample  mean  and  variance  of  the  feature  c  in 
class  A,,  respectively. 

Then,  it  should  be  noticed  that  class  A,-  could  consist  of  several  different  texture 
classes  in  the  case  that  the  different  textures  share  the  same  fractal  scale  (the 
roughness  of  the  surface),  but  have  different  patterns.  This  means  that  sometimes, 
checking  the  fractal  scale  only  is  not  enough  to  distinguish  the  different  patterns  of 
texture.  Thus,  we  need  an  additional  classification  scheme  to  distinguish  these  even  in 
the  same  class  A,.  This  additional  classification  will  be  done  in  second  level  of 
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classification  structure.  Here,  the  number  of  first-level  classes  M  is  chosen 
considering  the  sample  means  and  variances  of  the  fractal  scale  c  of  the  sample 
texture  images.  For  example,  if  two  different  texture  patterns  have  close  enough 
sample  means  of  fractal  scale  each  other,  such  that  the  range  [ c  -  cc  ,c  +  oc  ]  of  one 
texture  is  overlapped  with  the  other,  these  two  texture  patterns  could  be  classified  to 
the  same  class  in  first-level.  Thus,  it  should  be  noticed  that  it  is  possible  for  a  texture 
to  belong  to  several  different  classes  simultaneously,  depending  its  value  of  c. 

5.5.2.  The  Second-Level  of  Classification 

This  second-level  of  classification  is  an  optional  procedure.  If  one  fractal  scale 

feature  is  enough  to  classify  the  texture  in  first-level,  in  other  words,  if  the  class  At  in 

first-level  contains  only  one  member  in  it,  our  classification  could  be  terminated  in 

first-level.  On  the  other  hand,  if  the  class  has  several  members  in  it,  the  classification 

is  continued  to  the  second-level.  In  second-level,  the  textures  which  were  already 

classified  to  the  same  class  in  first-level  are  split  to  the  different  sub-classes,  based  on 

the  values  of  pattern  features  coi ,  CD 2  in  the  second-order  fractional  differencing 

periodic  function  (5.3.1.21).  Thus,  the  estimation  of  those  parameters  can  be 

completed  by  a  Maximum  Likelihood  Estimation  (MLE)  technique  discussed  in 

~  ( k )  -  ( k )  ~  ( k ) 

section  5.4.2,  and  these  extracted  features  denoted  by  Q  =  {(0\  ,o>2  }  are  used 

for  further  classification  by  measuring  another  weighted  distance  between  these 
features  and  the  mean  feature  of  each  of  the  N  subclasses  in  a  particular  class  At. 
Similarly  to  the  first-level, 

k*  =  minimum  (!(&  \k),  k  =  ],...,N  (5. 5. 2.1) 

k 


where 
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d(tik\k)=  X  f-j-f  f-]-  {5.522) 

/=“ l-«2  Xta/^) 

M 

and  /  and  [Gy  p  ;  correspond  to  the  sample  mean  and  variance  of  subclass  (k) 
features,  respectively.  Here,  it  should  be  noticed  that  since  we  have  at  most  several 
subclasses  from  a  first-level  class,  we  need  to  compare  only  a  small  number  of 
subclasses  to  complete  the  classification,  instead  of  checking  the  feature  distance  of 
whole  other  texture  classes.  Thus,  we  can  save  the  total  processing  time  of 
classification  by  using  this  multi-level  structure.  Therefore,  from  this  multi-level 
classification  method,  we  could  get  the  accurate  classification  and  save  the  processing 
time  at  the  same  time. 

5.6.  Experimental  Results 

For  these  experiments,  nine  different  classes  of  texture  were  taken  from 
Brodatz’s  standard  texture  album  [28]  for  the  training  set.  These  are,  namely, 
grass[D9],  tree  bark[D12],  straw[D15],  herringbone  weave[D17],  woolen  cloth[D19], 
calf  leather[D24],  beach  sand[D29],  water[D37],  and  raffia[D84],  Figure  5.4  shows 
the  256x256  original  texture  images  of  these. 

For  the  actual  training,  sixteen  64x64  sized  sample  image  data  were  taken  for 
each  different  texture  pattern,  and  the  sample  mean  and  variance  of  parameters,  c,  G>i , 
and  (jl>2  were  obtained  for  each  texture  class,  based  on  the  first  and  second-order 
fractional  differencing  models  [Table  5.1].  As  we  can  see  from  Table  5.1,  fractal 
scale  c  itself  is  not  enough  to  classify  the  different  textures,  because  some  of  textures 
have  similar  values  of  c,  even  though  they  are  different  texture  patterns.  This  is  the 
ieason  why  the  classification  will  be  completed  in  second  level  by  considering 


Figure  5.4:  256x256  original  texture  images  of  training  set:  (a)  grass  (b)  tree  bark 
(c)  straw  (d)  herringbone  weave  (e)  woolen  cloth  (0  calf  leather  (g)  beach  sand  (h) 
water  (i)  raffia 
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additional  parameter  values  of  <*>!,  a>2.  As  mentioned  before,  since  the  estimation  of  c 
can  be  done  by  simple  Least-square  estimation  and  the  estimation  of  CO] ,  CO2  can  be 
done  by  Maximum  Likelihood  estimation,  the  two-level  hierarchical  classification 
structure  induces  the  reasonably  reduced  processing  time  preserving  the  accuracy  of 
the  classification.  Based  on  these  sample  mean  and  variance  values,  the  number  of 
classes  for  the  first  level  of  classification,  M  =  5,  were  determined.  The  following  table 
[Table  5.2]  shows  the  final  classes  for  the  first  level  of  classification,  and  the 
corresponding  sample  mean  and  variance  f  each  class.  Notice  that  the  herringbon 
weave  texture  belongs  to  the  class  2  and  3,  because  of  its  high  value  of  variance. 
Then,  for  the  second  level  of  classification,  subclasses  can  be  taken  as  the  members  of 
each  class  based  on  the  different  values  of  0)1  and  0)2. 

5.6.1.  2-D  Rotated  Texture  Case 

In  this  experiment,  the  test  input  images  were  taken  from  the  2-D  raffia  textures 
rotated  by  various  angle  0s  (Figure  5.5).  Then,  each  64x64  texture  was  classified  by 
the  proposed  multi-level  classification  scheme.  For  the  first  level,  the  fractal  scale 
parameter  c  was  extracted  based  on  the  first-order  Fractional  Differencing  model 
(5.3.1),  and  the  parameters,  and  o>2,  were  extracted  from  the  second-order 
Factional  Differencing  periodic  model  (5.3.1.11).  Actual  classification  of  the  test 
images  was  done  in  each  level  by  comparing  weighted  distances  (5.5. 1.1 -2,  5.5.2. 1-2) 
between  the  extracted  features  and  the  data  base.  The  classification  results  are 
presented  in  Table  5.3.  Table  5.3  shows  the  parameter  values  extracted  from  each 
rotated  texture  pattern  and  the  demostrates  the  perfect  result  of  classification  based  on 


these  values. 


127 


Table  5.1:  The  sample  mean  and  variance  of  parameter  c,  (Oj,  02:  16  64x64 
sample  image  data  are  taken  for  each  different  texture  classes,  and  the  parameter 
values  are  extracted  from  the  first  and  second-order  fractional  differencing  models 
(5.3. 1-2). 


Textures 


grass 
tree  bark 
straw 

herringbone  weave 
woolen  cloth 
calf  leather 
beach  sand 
water 
raffia 


1.003  0. 

0.809  0. 

1.064  0. 

1.195  0. 

1.074  0. 

1.547 


0. 

0. 

0.062 


1.175 

0.665 

0.083 

1.042 


0.114  1.175 

0.095  0.793 


0.114 

0.107 

0.064 

0.153 


0.935 

0.571 

0.972 

0.988 


1 


.082 

.324 

.070 

0.167 

0.098 

0.122 

0.129 

0.132 

0.165 


Table  5.2:  Database  of  the  first  level  of  classification,  c,  and  a,2  are  the  sample 
mean  and  the  variance  of  class  /,  respectively. 


Class 


Textures 


woolen  cloth 

straw,  herringbon  weave 

herringbon  weave,  calf  leather,  water 

grass,  beach  sand 

tree  bark,  raffia 


0.055 
1.202  0.045 

1.539  0.067 


Table  5.3:  Classification  results  from  the  2-D  rotated  texture  images.  (Result 
indicates  the  result  class  after  applying  2-level  classification  method.) 


Angles 

c 

CO, 

.Ji>2 

Result 

■mm m 

■mw 

mwmm 

EMI 

1.517 

1.144 

1.102 

raffia 

Kfl 

1.535 

1.138 

1.119 

raffia 

■&< 

1.537 

1.142 

1.120 

raffia 

100 

1.532 

1.139 

1.118 

raffia 

120 

1.529 

1.138 

1.120 

raffia 

140 

1.527 

1.135 

1.097 

raffia 

160 

1.533 

1.140 

1.113 

raffia 

180 

1.525 

1.133 

1.099 

raffia 

5.6.2.  Orthographically  Projected  Texture  Case 

Six  64x64  test  input  images  were  taken  in  this  experiment  from  the  herringbon 
weave  textures  projected  orthographically  from  the  various  tilted  and  slanted  texture 
surface  (Figure  5.6).  Then,  each  texture  was  classified  by  the  proposed  multi-level 
classification  scheme.  Like  in  previous  experiment,  for  the  first  level,  the  fractal  scale 
parameter  c  was  extracted  based  on  the  first-order  Fractional  Differencing  model 
(5.3.1),  and  the  parameters,  (Dj  and  0)2,  were  extracted  from  the  second-order 
Factional  Differencing  periodic  model  (5. 3. 1.2).  The  parameter  values  extracted  from 
each  projected  texture  pattern  and  the  classification  results  from  this  experiment  are 
presented  in  Table  5.4.  Here,  notice  that  the  first  test  texture  was  assigned  to  Class  3 
for  the  first  level,  and  misclassified  to  the  calf  Ieathre  texture,  because  these  two 
texture  patterns  share  the  close  values  of  parameters. 
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Table  5.4:  Classification  results  from  the  orthographically  projected  herringbone 
weave  texture  images.  (Result  indicates  the  result  class  after  applying  2-level 
classification  method.) 


Angles 

'"x  =1T,  a  =T5^ 
x  =  0°,  a  =  30° 
x  =  0°,  o  =  45° 
x  =  45°,  o  =  15° 
x  =  45°,  a  =  30° 
x  =  45°,  a  =  45° 


c  I  CQi  I  to?  [  Result 

0.997  1.189  1 .0/7  calf  leather 

1.012  1.223  1.102  herringbone  weave 

1.055  1.218  1.096  herringbone  weave 

0.980  1.197  1.113  herringbone  weave 

1.023  1.209  1.100  herringbone  weave 

1.046  1.207  1.108  herringbone  weave 


mm  mmm 


i> 


5.6.3.  Rotated  And  Projected  Texture  Case 


In  this  experiment,  six  64x64  test  input  images  were  taken  from  the  straw 
textures  rotated  and  projected  orthographically  from  the  various  tilted  and  slanted 
texture  surface  (Figure  5.7).  Like  in  previous  experiments,  for  the  first  level,  the 
fractal  scale  parameter  c  was  extracted  based  on  the  first-order  Fractional  Differencing 
model  (5.3.1),  and  the  parameters,  £0i  and  CO2,  were  extracted  for  the  second  level  of 
classification  from  the  second-order  Factional  Differencing  periodic  model  (5.3.1.17). 
The  classification  results  from  this  experiment  are  presented  in  Table  5.5.  Table  5.5 
shows  the  parameter  values  extracted  from  each  rotated  and  projected  texture  pattern 
and  demonstrates  the  perfect  result  of  classification  based  on  these  values. 

5.7.  Conclusions 

A  new  multi-level  classification  technique  was  proposed  to  classify  the  3-D 
rotated  texture  surface.  Unlike  most  classification  schemes  which  have  been  suggested 
to  date,  this  classification  method  can  handle  arbitrary  3-D  rotated  samples  of  textures, 
i.e.,  the  accuracy  of  classification  is  not  affected  by  the  3-D  rotation  of  the  test  texture. 
The  proposed  multi-level  classification  scheme  consists  of  two  levels  of  classification 
procedure.  In  the  first  level  of  classification,  a  3-D  rotational  invariant  feature  c 
(fractal  scale)  in  the  first-order  Fractional  Differencing  model  was  extracted,  and 
based  on  this  value,  the  test  data  image  was  classified  to  a  certain  class.  In  the  second 
level,  each  class  was  divided  to  the  final  desired  subclasses  based  on  two  other  texture 
pattern  features,  (Oj  and  (1)2,  which  were  extracted  from  the  second-order  Fractional 
Differencing  periodic  model.  Then  the  input  texture  image  was  classified  further  in 
detail  with  these  two  pattern  features.  This  multi-level  classification  has  several 
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Table  5.5:  Classification  results  from  the  rotated  and  orthographically  projected 
straw  texture  images.  (Result  indicates  the  result  class  after  applying  2-level 
classification  method.) 


A 


advantages  over  the  conventional  one-level  classification  methods.  First,  since  this 
algorithm  uses  a  simple  LS  estimation  to  get  the  fractal  scale  c  in  first  level,  the 
processing  time  of  the  classification  is  reduced  comparing  to  the  methods  that  rely  on 
only  ML  estimation.  Second,  from  the  additional  texture  patterns  C0i ,  0)2  in  second 
level,  the  more  accurate  classification  can  be  achieved  compared  to  the  classification 
methods  which  have  only  one  feature  parameter  c.  As  a  result  of  a  series  of 
experiments  involving  the  differently  oriented  samples  of  natural  textures,  it  is 
concluded  that  these  combined  features  make  possible  for  this  multi-level 
classification  method  to  have  a  strong  class  separability  power  for  arbitrary  oriented 
3-D  texture  patterns. 
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CHAPTER  6 
CONCLUSION  AND 
FUTURE  RESEARCH 


6.1.  Conclusions 

The  modeling  and  analysis  of  texture  pattern  of  a  three  dimensional  surface  have 
been  investigated.  Since  3-D  natural  scene  image  contains  both  radiance  and  texture 
information,  3-D  texture  analysis  can  not  be  done  by  either  shading  analysis  or  texture 
analysis  only.  Also,  the  distortion  of  the  regular  texture  pattern  due  to  the  3-D  surface 
orientation  makes  it  difficult  to  analyze  3-D  texture  by  the  conventional  stationary 
random  field  texture  model.  Thus,  in  this  thesis,  3-D  surface  model  which  can  handle 
both  shading  and  texture  information  and  3-D  texture  analyses  which  are  not  affected 
by  3-D  surface  orientation  have  been  emphasized. 

The  contribution  of  the  research  can  be  summarized  as  follows.  First,  the 
orthographically  projected  fractional  differencing  model  was  introduced  to  represent 
the  orthographically  projected  texture  pattern  due  to  the  3-D  surface  orientation.  This 
model  has  an  ability  to  synthesize  a  long-correlated  or  a  short-correlated  random 
texture  with  various  values  of  fractal  scale  parameters  c  and  d,  and  the  roughness  and 
the  distorted  texture  pattern  of  the  3-D  surface  simultaneously. 

Second,  the  estimation  scheme  of  the  parameters  of  the  orthographically 
projected  fractional  differencing  model  was  developed  based  on  a  hybrid  method  of 
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the  least  mean  square  and  maximum  likelihood  estimations.  For  the  estimation  of  the 
fractal  scaling  parameter,  this  estimation  scheme  is  much  simpler  than  the  fractional 
Brownian  process  based  model’s,  because  of  its  discrete  process.  Rotational  and 
scaling  invariant  properties  of  the  parameters  were  presented.  Thus,  even  after 
performing  several  linear  coordinate  transformations,  all  Gaussian  assumptions  of 
random  noise  were  preserved. 

Third,  a  composite  model  of  the  Shape-from-shading  and  the  Shape-from-texture 
has  been  proposed  to  represent  a  natural  scene  which  contains  both  shading  and 
texture  information,  considering  the  scene  as  the  superposition  of  a  deterministic 
function  and  a  random  texture.  Surface  orientations  and  the  illumination  direction 
were  extracted  directly  from  a  single  natural  scene  image  without  any  pre-processing. 
This  avoids  the  possible  error  cumulation  from  each  process  step.  Also,  comparing 
with  the  conventional  Shape-from-shading  techniques,  this  composite  model  gives 
unique  and  more  reliable  solutions  for  the  surface  orientation  parameter  values, 
because  of  the  additional  constraint  from  the  texture  function  part. 

Fourth,  a  multi-level  3-D  rotational  invariant  classification  scheme  has  been 
developed,  based  on  the  first  and  second-order  fractional  differencing  models.  Multi¬ 
level  structure  of  this  classification  consists  of  two  hierarchical  levels.  In  the  first 
level,  the  fractal  scale  which  is  known  to  be  a  rotational  and  scaling  invariant 
parameter  was  extracted  from  the  first-order  fractional  differencing  model,  and  with 
this  value,  preliminary  classification  has  been  done.  In  the  second  level,  the  more 
detailed  classification  has  been  achieved  with  additional  texture  pattern  parameters 
from  the  second-order  fractional  differencing  model.  This  hierarchical  structure  could 
reduce  the  total  processing  time  from  the  simple  least  mean  square  estimation  scheme 
for  the  first  level,  and  at  the  same  time,  it  could  preserve  the  accuracy  of  the 
classification  from  the  second  level  of  classification. 
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6.2.  Suggestions  For  Future  Research 

In  addition  to  the  pattern  analyses  of  the  texture  on  3-D  surface,  which  are 
developed  in  this  thesis,  several  related  works  can  be  suggested  as  follows. 

6.2.1.  Neural  Network  Analysis  For  The  Global  Optimization  Problem 

In  the  maximum  likelihood  estimation  procedure,  the  optimization  routine 
maximizes  the  likelihood  function.  Since  our  likelihood  function  is  a  complicated 
non-linear  function  and  the  existing  approximation  techniques  for  getting  the 
maximum  value  from  the  non-linear  function  require  a  great  many  iterations,  the 
computational  time  is  enormous.  Recently,  based  on  the  parallel  processing  structure 
of  the  computer,  some  researchers  claim  that  neural  network  analysis  is  one  of  the 
promising  techniques  for  reducing  the  computational  time  for  optimization  [163]. 
Therefore,  investigating  the  applicability  of  neural  network  analysis  to  the 
maximization  of  the  likelihood  function  will  be  valuable. 

6.2.2.  Modified  Scheme  Of  The  Shape  From  Shading  And  Texture 
Method  With  Weighting  Factor 

The  proposed  composite  model  of  the  Shape-ffom-shading  and  the  Shape-from- 
texture  does  not  have  any  weighting  factor  between  both  3-D  analysis  techniques. 
However,  in  situation  that  the  radiance  information  is  more  dominant  than  the  texture 
information,  or  in  the  reverse  situation,  the  estimated  surface  orientations  from  the 
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Shape-from-shading  or  Shape-from-texture  technique  only  can  be  more  reliable. 
Therefore,  if  we  have  a  proper  criterion  to  select  either  technique  or  the  proposed 
composite  model  to  represent  the  part  of  the  surface  image  data,  we  may  improve  the 
results.  One  of  the  possible  way  to  do  this  will  be  given  by  checking  the  fractal  scale 
of  the  surface.  In  [176],  Pentland  described  the  fractal  scaling  factor  as  a  measure  of 
the  roughness  of  surface.  That  is,  if  the  fractal  scale  c=0,  it  represents  the  smooth 
surface  with  gentle  random  undulations.  In  contrast,  the  surface  with  fractal  scale  c>0 
are  not  perceived  as  smooth,  but  rather  as  being  rough  or  three-dimensionally 
textured.  Thus,  by  checking  the  estimated  value  of  c,  it  can  be  possible  to  make 
decision  which  technique  should  be  preferred. 

6.2.3.  High  Resolution  Analysis  For  The  Texture  Pattern 

Our  proposed  texture  model  in  chapter  4  is  based  on  local  patch  analysis.  In  the 
simulation  part  of  this  thesis,  we  used  a  16  x  16  sized  patch  taken  from  the  center  of  a 
32  x  32  sized  patch  to  build  the  whole  image.  There  is  a  minimum  requirement  on  the 
size  of  each  patch  for  being  able  to  extract  meaningful  statistical  information  from  the 
patch.  This  is  the  bottleneck  for  the  high  resolution  analysis.  Since  the  elementary 
unit  for  surface  orientation  is  a  single  patch  (we  assume  each  patch  to  be  a  plane  and 
consequently  the  values  for  the  surface  orientation  parameters  a,  x  are  constants  over 
the  patch),  the  end  result  is  a  low  resolution  analysis.  Therefore,  we  need  to  develop  a 
high  resolution  scheme  based  on  a  single  pixel  instead  of  a  complete  patch.  Therefore, 
if  we  can  develop  a  pixel  based  model,  the  synthesis  of  3-D  texture  can  be  done  by 
just  adopting  the  surface  orientation  parameter  values  recursively  for  each  pixel. 
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6.2.4.  Robust  Estimation  For  The  Contaminated  Gaussian  Noise  Case 


The  estimation  procedure  which  is  proposed  in  this  thesis  is  based  on  the 
assumption  that  the  input  noise  has  a  pure  Gaussian  distribution.  However,  in  the  real 
world,  pure  Gaussian  noise  does  not  exist.  Therefore,  even  if  we  have  noise  which  is 
made  of  Gaussian  noise  mixed  with  a  small  portion  of  other  distributed  noise,  our 
estimation  will  not  give  the  right  estimated  values  for  the  parameters.  Therefore,  it  is 
favorable  to  develop  a  robust  estimation  scheme  for  this  situation. 

6.2.5.  Boundary  Detection  Of  The  Mixed  3-D  Textures  With  Fractal  Scale 
And  Knowledge-Based  Post-Cleaning  Process 

The  natural  scene  image  does  not  contain  a  pure  2-D  texture  pattern  which  can 
be  represented  by  the  isotropically  distributed  random  texture  model,  but  a  projected 
pattern  on  the  3-D  surface.  For  example,  in  Figure  2.2  (tree  image),  the  tree  bark 
pattern  around  the  boundary  between  the  tree  and  the  lawn  ground  looks  more  dense 
than  the  pattern  in  the  middle  of  the  tree.  This  makes  it  difficult  to  apply  the 
conventional  2-D  texture  model  to  get  the  texture  boundary  the  different  textures. 
However,  this  3-D  texture  boundary  can  be  detected  using  fractal  scaling  parameter  in 
fractional  differencing  model  which  has  been  proposed  through  this  thesis  because  the 
fractal  scale  parameter  is  known  to  be  a  rotational  and  scaling  invariant  parameter  and 
its  value  is  considered  as  a  measure  of  the  roughness  of  surface.  Thus,  we  can  have  the 
same  value  of  the  fractal  scale  within  same  3-D  texture  pattern,  and  the  higher  value 
of  the  fractal  scale  around  the  texture  boundary  [46, 173],  Then,  the  fractal  scale 
image  can  be  generated  by  applying  a  proper  size  of  window  to  the  original  image, 
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and  the  texture  boundary  can  be  obtained  by  thresholding  the  values  of  fractal  scale  in 
the  fractal  scale  image.  Also,  this  texture  boundary  image  can  be  improved  by 
deleting  false  boundaries,  correcting  wrong  boundary  directions,  and  linking  the  lost 
boundaries.  This  post-processing  can  be  done  by  applying  a  knowledge-based 
reasoning  process  to  the  obtained  boundary  image.  Here,  a  set  of  rules  for 
knowledge-based  algorithm  can  be  constructed  from  the  prior  knowledge  that  the 
texture  boundaries  are  continuous  and  closed.  This  boundary  detection  technique  is 
being  developed  in  detail. 
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Appendix  A 

Proof  of  Theorem  2.1 

Intensity  function  x(l]yl2)  can  be  represented  by  the  following  equation. 

x(lx, l2)  =  X\i(L-N)  (A.l) 

where,  X|i:  albedo  (constant  for  the  rcflecuiess) 

L:  the  illumination  direction  vector 
N:  the  surface  normal  vector 

Ihus,  when  (/ \,l2)  defines  the  image  plane  from  the  viewing  direction,  the  first 
derivative  of  image  intensity  in  the  direction  {dl\,dl 2),  dx{l\,l2)  is 

dx  =  \|i(dL-N  +  L-dN)  (A.2) 

Assuming  L  is  constant, 

dx  =  }qt(L'dN)  (A. 3) 

Notice  that  d N,  the  change  in  N,  is  perpendicular  to  N  as  it  lies  in  the  tangent  plane  to 
N,  and  isotropically  distributed  as  same  as  N  is.  Thus 

£(IdN)  =  0  (A.4) 

and 

E(dx)  =  (dN))  (A.5) 

+LhdNh)  (A.6) 

where  (dN^  ,dNi2,dN ;3)  are  the  average  values  of  change  in  the  surface  normal  in 
image  direction  (d/ 1  ,dl2 )  and  L=(L/(  ,/_/2,L/3)  is  the  illumination  direction.  That  is, 


159 


Since  dN^  =  0, 

E(dx)  ~  X|i(L/id!/V/1  +  L[2dNi2)  (A. 8) 

Then,  introducing  dr,  which  may  be  thought  of  as  the  expected  magnitude  of  dN,  that 
is,  as  £(|dN|), 

dNh=l\dr,  dNtl  =  l2dr  (A-9) 


where  /l2+/22=  1  and  lxll2=dfiiJMll=dxiJdxll,  defining  (dxlrdxh)  to  be  the 
differential  step  in  the  image  along  which  dx  was  measured. 

Thus,  defining  Lt=X\iLixdr  and  Lt=X[iLl2dr,  we  can  have  the  following  linear 
regression  model,  from  (A.7),(A.8)  and  (A.9). 


dl  ii  dl 2i 

dl\2  dl 22 

dl  \n  dl  In 


(A.  10) 


where  dxx  is  the  average  value  of  dx  over  the  i-th  patch  in  direction  (dl  u,dl2i)- 
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Appendix  B 

Proof  of  Theorem  2.2 

If  we  denote  the  matrix  of  directions  (dlu,dl2i)  by  P,  and  let  pT  indicate  the 
transpose  of  P,  then  the  solution  (A.  10)  in  Appendix  A  is  the  following  least-square 
estimator. 

dx\ 

nr  ^ 

=  #pr'pr  . 

dXn 

Here,  notice  that  the  above  solution  is  a  normalized  one.  Thus,  in  order  to  get  the 
actual  values  of  ( L ,L/2,L/3),  we  need  to  calculate  the  value  of  X^dr  properly. 

E(dK2)  =  E(lV(LlldNli  +LlldNll  +LlidNli)2)  (B.2) 

=\2\i2{Ltl2dNlx2  +Lh2dNlt2  +L,2dNl\2Ll{LlldNlldNh 

+  2LliLljdNlidN,i+2L,lL,3dNl2dNlj)  (B.3) 

where  ( dN ^  ,dN[2,dN[2)  are  the  average  values  of  (dN^  ,dNi2,dNi3). 

Since  Var {dN^)  =  Var(d7vf/2)  =  dr 2  for  sphere  model,  dN^2  =  l2dr2jcdr2,  dN\2  = 
l22dr2+  dr 2,  and  dN )3  =  dr2, 

E{dx2)  -  l2\i2(L[2(ll2dr2+dr2)+Ll2(l22dr2+dr2) 

+L,2dr1+2L,xLl2lll2dr2)  (B.4) 


(B.l) 
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=  XV2((/i^l+/2^/2)2^2+(^12+^22+^/,2)^2)  (B.5) 

If  the  illumination  direction  vector  L  is  counted  as  a  unit  vector,  that  is,  Li2  +  Li2  + 

^/22  =  1, 

E{dx2)  =  \2\i2{{lxLlx+l2Lh)2dr2+dr2)  (B.6) 

Also,  from  (2.2.1.1.7),(2.2.1.1.8),  and  (2.2. 1.1. 9), 

E{dx)2  -{X^dNi^+dNi^))2 

-(kiM^+^LiJdr)2 

-XVC/^+Z^/dr2  (B.7) 

Thus, 

E(dx2)-E(dx)2 -X2|i2dr2  4  k2  (B.8) 

Therefore,  from  the  definitions  of  ,  L/2,  and  the  equation  L/]  2  +  L{2  +  L[2  =  1,  the 
illumination  direction  can  be  calculated  as  follows. 


And,  from  the  relation  between  the  x-y-z  coordinate  system  and  the  angular  coordinate 
system,  we  can  represent  the  illumination  direction  with  its  tilt  and  slant  angles,  x^, 
<?L.  Thus, 

i  4  , 

*L  —  tan  (-x“),  O/,  =cos  lL /3 
Llx 


(B.10) 
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Appendix  C 

Proof  of  Theorem  2.3 

Let  be  the  phase  angle  between  the  viewing  direction  V  and  the  illumination 
direction  L,  y  be  the  incident  angle  between  the  illumination  direction  and  the  surface 
normal  direction,  and  T|  be  the  emittance  angle  between  the  viewing  direction  and  the 
surface  normal  direction.  Then 

V-L  =  cosa^,  L-N  =  cosy,  and  V-N  =  cos,n  (C.l) 

Intensity  function  x(l \,l2)  can  be  represented  by  the  following  equation. 

x(l \,lj)  —  ^M-(L'N)  =  Xpcosy,  when  cosy>0.  (C.2) 

Then,  for  each  of  surface  patches,  to  determine  the  expected  value  of  image  intensity, 
we  need  to  average  x (/ 1  ,/2)»  and, to  determine  the  expected  value  of  image  intensity 
squared,  we  need  to  average  jt2(/ i  ,l2)  over  the  image  of  a  hemisphere. 

Notice  that  since  not  all  of  the  hemisphere  is  illuminated,  we  should  integrate 
only  over  the  illuminated  region,  to  avoid  areas  where  cosy  is  negative.  A  suitable 
spherical  coordinate  system  can  be  erected  with  the  equator  in  the  plane  containing  L 
and  V,  with  the  pole  at  right  angles  to  this  plane,  namely  at  O,  defined  to  be  a  unit 
vector  orthogonal  to  the  plane  containing  L  and  V.  Let  latitude  be  x  measured  from 
the  equator,  while  longitude  is  c,  measured  along  the  equator  from  the  point  V  toward 
L.  Therefore,  the  point  L  has  longitude  and  zero  latitude.  Simple  spherical 
trigonometry  using  a  triangle  with  comers  V,  N,  and  O,  yields 

COST|  =  COST  COSO  (C.3) 


and  from  a  similar  triangle  with  comers  L,  N,  and  O,  we  can  obtain 
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cosy  =  cost  cos  (o  -  a/J  (C.4) 

The  illuminated  half  of  the  sphere  runs  from  CT  =  —k/2+Ol  to  ct  =  +7C/2-Kj^,  while  the 
visible  half  goes  from  a  =  -k/2  to  a  =  +n/2  and  the  infinitesimal  element  of  area  is 
cosxdadt. 

The  integrals  that  we  are  interested  in  are  of  the  form 
n/2  n/2 

J  j  /( g,t)  cost)  cost  dadt  (C.5) 

-TV 2  -TU2+oL 

where  the  cost]  term  compensates  for  the  foreshortening  due  to  the  projection  of  the 
spherical  surface  into  the  image.  We  need  this  factor  because  we  are  using  a 
coordinate  system  on  the  hemisphere,  but  are  seeking  an  average  over  the  image  of 
the  hemisphere.  If  we  include  self-shadowed  areas  in  the  computation  of  the  average, 
we  must  divide  the  integral  by  the  whole  area,  %,  of  the  disc  that  is  the  projection  of 
the  hemisphere,  on  the  other  hand,  if  we  do  not  include  self-shadowed  areas,  we 
divide  by  the  area  ( n/2)(l+coso i)  of  the  projection  of  the  illuminated  part  of  the 
hemisphere. 

Thus,  from  (C.2)  and  (C.4), 

E  (x)  =  A.|i£  (cosy) 

=  X.p£  (cost  cos  (a  -  c^))  (C.6) 

Here,  to  get  the  value  of  E  (cosy),  we  need  to  evaluate  the  integral 

71/2  71/2 

B  i  =  j  j  (COST]  cosy)  cost  dodx  (C.7) 

-tv2  -71/2+0, 
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n/2  ti/2 

=  J  J  (cos x  caso)(ty;.vx  cos  (o-oLj)  cost  dodx  (C.8) 

-ti/2  -ti/2 fO/, 


71/2  71/2 

=  |  cos3xdx  j  cos  (o-Ol)  coso  do 

-n/2  -tv'2+ol 


(C.9) 


=  ~(sinoi+(n-oL)cosoL) 


(C.10) 


Therefore, 


E  (x )  —  ~~(sin oL+(%-oL )cos ol), 
3  k 


(C.ll) 


E  (x)  =  - - -(sinoL+(n-oL)cosoL), 

3k(  1+mvci/J 


(C.  12) 


depending  on  whether  we  average  over  all  image  regions,  including  self-shadowed 
parts,  or  not. 

Similarly,  to  get  the  value  of  E  ( x 2),  we  need  to  evaluate  the  integral 


n/2  n/2 


/#2=  j  j  cost|  cos^t}/  cost  dodx 

-n/2  -R/2fO/. 


(C.13) 


7t/2  n/2 


J  J  (cos  t  coso)(cosx  cos  (o-oL))2cosx  dodx  (C.14) 

n/2  ;i/2-tCT,. 


j  cos4xdx  j  cos2 (o-Ol)coso  do 
n/2  n/2+o, 


(C.  15) 


-( l+eoscT/ X 


(C.  16) 


Therefore, 
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£(*2)=-^~ ^ -d+cosaL)2,  (C.17) 

O 

or 

£(x2)  =  -^-(1+cosol),  (C.18) 

depending  on  whether  we  average  over  all  image  regions,  including  self-shadowed 
parts  or  not. 
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